Method for predicting propagation of sound in an ocean

ABSTRACT

Propagation of sound in an ocean is accurately predicted at all ranges from a source by a computational method that uses a newly created Decomposition Method to make the predictions at long range from the source. The predictions are done in a limited region of the ocean without the need to predict propagation of sound outside the region. The vertical boundaries of the region are parallel to a given direction of propagation of sound from the source. These boundaries are transparent, and sound penetrates them without reflection. The Decomposition Method is based on a special splitting of the sound field into primary and secondary components. The primary components are related to the sound field that would exist in the region if the vertical boundaries of the region were acoustically hard. They are found by integrating a system of coupled parabolic equations that are derived from the coupled mode equations for the sound field in the region. The primary components contain spurious reverberations that originate at the vertical boundaries of the region. The secondary components identify and cancel the spurious reverberations in the primary components. The secondary components are found by integrating standard parabolic equations subject to transparent boundary conditions at the vertical boundaries of the region. These boundary conditions also include source terms that are obtained from the primary components. When the primary and secondary components are summed, the spurious reverberations are cancelled in the sum. The actual sound field in the region is found from this sum. A modification of the Decomposition Method incorporates an extrapolation procedure to ensure accuracy of the predictions at long range from the source.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of PPA Ser. No. 60/478,573, filed 2003 Jun. 12 by the present inventor.

FEDERALLY SPONSORED RESEARCH

Not applicable

SEQUENCE LISTING OR PROGRAM

Not applicable

BACKGROUND OF THE INVENTION

1. Field of the Invention

The invention pertains to a method for predicting propagation of sound in an ocean.

2. Description of the Related Art

For the purpose of predicting propagation of sound in an ocean, the ocean can be modelled as a three-dimensional acoustic waveguide that is spatially irregular. A partial differential wave equation governs propagation of sound in the waveguide. A complex-valued distribution of pressure represents the sound field in this equation. However, the sound field also can be represented in terms of the local modes of the waveguide. This representation is useful because it reduces the dimensionality of the model. The equations that govern this representation are obtained by taking moments of the wave equation with respect to the local modes. These moments produce a system of coupled mode equations in two horizontal dimensions. The solution of these coupled mode equations can be used to reconstruct the sound field at any depth in the waveguide from the local mode representation.

In FIG. 1 there is shown a circular region 10 in the horizontal plane with a source of sound at the center 12. An article by A. T. Abawi, W. A. Kuperman and M. D. Collins (“The coupled mode parabolic equation”, Journal of the Acoustical Society of America, Vol. 102, No. 1, pp. 233-238, July 1997) discusses a numerical method for solving coupled mode equations to predict propagation of sound over region 10. This method is an elaboration of an earlier numerical method that was developed by M. D. Collins ignoring coupling among the modes. Propagation of sound is predicted over all of region 10. The preferred direction of propagation is the radially outward direction from the center 12.

An article by M. D. Collins (“The adiabatic mode parabolic equation”, Journal of the Acoustical Society of America, Vol. 94, No. 4, pp. 2269-2278, October 1993) points out that it would be very efficient to make predictions just over an angular sector 14 of small angular width. The article notes that transparent boundary conditions would have to be imposed along the radial edges of sector 14 for this purpose, but it does not indicate how to formulate these conditions or how to adopt them in a solution method. The article by A. T. Abawi, W. A. Kuperman and M. D. Collins does not address these issues at all.

SUMMARY OF THE INVENTION

The invention permits propagation of sound in an ocean to be accurately predicted at all ranges from a source. At long ranges this prediction is done by a newly created Decomposition Method. The predictions are done for fully three-dimensional ocean environments.

The invention permits the predictions to be done in a limited region of the ocean without the need to predict propagation of sound outside the region. The vertical boundaries of the region are parallel to a given direction of propagation of sound from the source. These boundaries are transparent, and sound penetrates them without reflection.

The Decomposition Method is based on a special splitting of the sound field into primary and secondary components. The primary components are related to the sound field that would exist in the region if the vertical boundaries of the region were acoustically hard. They are found by integrating a system of coupled parabolic equations that are derived from the coupled mode equations for the sound field in the region. The primary components contain spurious reverberations that originate at the vertical boundaries of the region.

The secondary components identify and cancel the spurious reverberations in the primary components. The secondary components are found by integrating standard parabolic equations subject to transparent boundary conditions at the vertical boundaries of the region. These boundary conditions also include source terms that are obtained from the primary components. When the primary and secondary components are summed, the spurious reverberations are cancelled in the sun. The actual sound field in the region is found from this sum.

A modification of the Decomposition Method incorporates an extrapolation procedure to ensure that the predictions are accurate at long range from the source.

It is an object of the invention to provide an accurate method for predicting propagation of sound in an ocean.

Another object of the invention is to provide an accurate method for predicting propagation of sound from a source in a limited region of an ocean without the need to predict propagation of sound outside the region.

Yet another object of the invention is to provide an accurate method for predicting propagation of sound from a source in a limited region of an ocean without the need to predict propagation of sound outside the region when the length of the region is much greater than the width of the region.

Other objects, advantages and new features of the invention will be apparent from the following detailed description of the invention when considered in conjuction with the accompanying drawings.

DRAWINGS

FIG. 1 shows a circular region where prior art methods predict propagation of sound.

FIG. 2 shows a section of a waveguide that models the ocean environment.

FIG. 3 shows initial and boundary conditions in the horizontal plane of range and cross-range.

FIG. 4 shows a flowchart for the preferred embodiment of the method of the invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

1. Coupled mode equations in slab geometry and in spherical geometry.

The propagation of sound in an irregular three-dimensional waveguide can be described mathematically by a partial differential equation subject to certain conditions at the boundaries of the waveguide. (By “irregular” we mean that the medium in the waveguide is spatially inhomogeneous. We also assume that this medium can be modelled as a lossless fluid.) Let the variables x, y and z denote range, cross-range and depth in the waveguide, respectively. Let the variables ω, c and ρ denote source frequency, local sound speed and local fluid density in the waveguide, respectively. We let the variable p denote the complex acoustic pressure. This quantity is defined by the property that the physical acoustic pressure is the real part of the product p exp(−iωt), where t denotes time. The partial differential equation in three dimensions that governs p is $\begin{matrix} {{{\frac{\partial\quad}{\partial x}\left( {\frac{1}{\rho}\frac{\partial p}{\partial x}} \right)} + {\frac{\partial\quad}{\partial y}\left( {\frac{1}{\rho}\frac{\partial p}{\partial y}} \right)} + {\frac{\partial\quad}{\partial z}\left( {\frac{1}{\rho}\frac{\partial p}{\partial z}} \right)} + {\frac{1}{\rho}\frac{\omega^{2}}{c^{2}}p}} = 0.} & (1.1) \end{matrix}$

Referring now to FIG. 2, we assume that the pressure vanishes at a reflective top surface 20 of the waveguide, where z=0, and at a reflective bottom surface 22 of the waveguide, where z=−L: p=0 at z=0 and at z=−L.  (1.2) We also assume that there is a single interface 24 between the surface and the bottom where the local sound speed and the local fluid density can be discontinuous. The depth of this interface can vary with range and with cross-range. We describe the interface by the equation S(x,y,z)=H(x,y)+z=0,  (1.3) where 0<H(x, y)<L. The acoustic pressure must be continuous across this interface. The component of acoustic fluid velocity that is normal to the interface must also be continuous across the interface. We can express these continuity conditions as follows. $p\quad{and}\quad\frac{1}{\rho}\left( {{\frac{\partial H}{\partial x}\frac{\partial p}{\partial x}} + {\frac{\partial H}{\partial y}\frac{\partial p}{\partial y}} + \frac{\partial p}{\partial z}} \right)$  are continuous across z=−H(x, y).  (1.4)

It is possible to convert the acoustic equation (1.1) and the conditions (1.2) and (1.4) on it to an equivalent infinite system of coupled partial differential equations in two dimensions. This system is based on the local modes of the waveguide, which are characterized by the eigenfunctions ƒ_(n)(x, y, z) and by the squared wavenumbers ξ_(n) ²(X, y) that satisfy the differential equation ${{\frac{\partial\quad}{\partial z}\left( {\frac{1}{\rho}\frac{\partial f_{n}}{\partial z}} \right)} + {\frac{1}{\rho}\left( {\frac{\omega^{2}}{c^{2}} - \xi_{n}^{2}} \right)f_{n}}} = 0$  if −L<z<0 and z≠−H(x,y),  (1.5) subject to the conditions that

ƒ_(n)=0 at z=−L and at z=0,  (1.6) $f_{n}\quad{and}\quad\frac{1}{\rho}\frac{\partial f_{n}}{\partial z}$  are continuous across z=−H(x, y).  (1.7)

At each position (x,y) these equations describe a nonsingular selfadjoint characteristic value problem. The squared wavenumbers ξ_(n) ²(x, y) are real and countably infinite in number. They can be indexed to decrease monotonically to −∞ as the index n increases through the positive integers. The local modes can be computed with numerical methods that are well-known in the art.

At each position (x, y), σƒ_(n)(x, y 0)/σz=0 if and only if ƒ_(n)(x, y, z) vanishes identically for −L≦z≦0. Therefore, we can assume that σƒ_(n)(x, z, y, 0)/σz≠0 and that ƒ_(n) (x, y ,z) does not vanish identically. Furthermore, we can assume that ƒ_(n)(x, y, z) is real-valued. Since the density ρ is positive, the eigenfunctions are orthogonal in the following sense: $\begin{matrix} {{\int_{- L}^{0}{\frac{1}{\rho}f_{m}f_{n}{\mathbb{d}z}}} = \left\{ \begin{matrix} {{{{d_{n}\left( {x,y} \right)} > {0\quad{if}\quad m}} = n},} \\ {{0\quad{if}\quad m} \neq {n.}} \end{matrix} \right.} & (1.8) \end{matrix}$

We also assume that the eigenfunctions are normalized so that the signs of the partial derivatives σƒ_(n)(x, y, 0)/θz and the values of the integrals d_(n)(x, y) are independent of the position (x, y).

The functions p, σp/σx and σp/σy have the following eigenfunction expansions: $\begin{matrix} {{p = {{\sum\limits_{n = 1}^{\infty}{a_{n}{f_{n}/d_{n}}\quad{where}\quad a_{n}}} = {\int_{- L}^{0}{\frac{1}{\rho}{pf}_{n}{\mathbb{d}z}}}}},} & (1.9) \\ {{\frac{\partial p}{\partial x} = {{\sum\limits_{n = 1}^{\infty}{b_{n}{f_{n}/d_{n}}\quad{where}\quad b_{n}}} = {\int_{- L}^{0}{\frac{1}{\rho}\frac{\partial p}{\partial x}f_{n}{\mathbb{d}z}}}}},} & (1.10) \\ {\frac{\partial p}{\partial y} = {{\sum\limits_{n = 1}^{\infty}{c_{n}{f_{n}/d_{n}}\quad{where}\quad c_{n}}} = {\int_{- L}^{0}{\frac{1}{\rho}\frac{\partial p}{\partial y}f_{n}{{\mathbb{d}z}.}}}}} & (1.11) \end{matrix}$ The coefficients a_(n), b_(n) and c_(n) in these expansions are functions of the position (x,y). They are governed by the following system of coupled partial differential equations. $\begin{matrix} {{{\frac{\partial b_{n}}{\partial x} + \frac{\partial c_{n}}{\partial y} + {\xi_{n}^{2}a_{n}}} = {\sum\limits_{m = 1}^{\infty}{\left( {{b_{m}B_{mn}} + {c_{m}C_{mn}}} \right)\frac{1}{d_{m}}}}},} & (1.12) \\ {{b_{n} = {\frac{\partial a_{n}}{\partial x} + {\sum\limits_{m = 1}^{\infty}{a_{m}B_{nm}\frac{1}{d_{m}}}}}},} & (1.13) \\ {c_{n} = {\frac{\partial a_{n}}{\partial y} + {\sum\limits_{m = 1}^{\infty}{a_{m}C_{nm}{\frac{1}{d_{m}}.}}}}} & (1.14) \end{matrix}$ These equations constitute the fundamental first-order formulation of coupled mode theory for an irregular three-dimensional waveguide.

The coupling coefficients B_(mn) and C_(mn) can be computed as follows. Let $\begin{matrix} {g_{n} = {\frac{1}{\rho}\frac{\partial f_{n}}{\partial z}}} & (1.15) \end{matrix}$ and define the following matrix elements: $\begin{matrix} {{D_{mn} = {\int_{- L}^{0}{\left\lbrack {{\frac{\partial\quad}{\partial x}\left( {\frac{1}{\rho}\frac{\omega^{2}}{c^{2}}} \right)f_{m}f_{n}} + {\frac{\partial\rho}{\partial x}g_{m}g_{n}}} \right\rbrack{\mathbb{d}z}}}},} & (1.16) \\ {{E_{mn} = {\int_{- L}^{0}{\frac{\partial\quad}{\partial x}\left( \frac{1}{\rho} \right)f_{m}f_{n}{\mathbb{d}z}}}},} & (1.17) \\ {{F_{mn} = {\int_{- L}^{0}{\left\lbrack {{\frac{\partial\quad}{\partial y}\left( {\frac{1}{\rho}\frac{\omega^{2}}{c^{2}}} \right)f_{m}f_{n}} + {\frac{\partial\rho}{\partial y}g_{m}g_{n}}} \right\rbrack{\mathbb{d}z}}}},} & (1.18) \\ {{G_{mn} = {\int_{- L}^{0}{\frac{\partial\quad}{\partial y}\left( \frac{1}{\rho} \right)f_{m}f_{n}{\mathbb{d}z}}}},} & (1.19) \\ {{J_{mn} = {{\left\lbrack {\frac{1}{\rho}\frac{\omega^{2}}{c^{2}}} \right\rbrack_{{- H} - 0}^{{- H} + 0}{f_{m}\left( {- H} \right)}{f_{n}\left( {- H} \right)}} + {\lbrack\rho\rbrack_{{- H} - 0}^{{- H} + 0}{g_{m}\left( {- H} \right)}{g_{n}\left( {- H} \right)}}}},} & (1.20) \\ {K_{mn} = {\left\lbrack \frac{1}{\rho} \right\rbrack_{{- H} - 0}^{{- H} + 0}{f_{m}\left( {- H} \right)}{{f_{n}\left( {- H} \right)}.}}} & (1.21) \end{matrix}$ If m≠n then $\begin{matrix} {{B_{mn} = {- {\left( {\xi_{m}^{2} - \xi_{n}^{2}} \right)^{- 1}\left\lbrack {D_{mn} - {E_{mn}\xi_{n}^{2}} + {\frac{\partial H}{\partial x}\left( {J_{mn} - {K_{mn}\xi_{n}^{2}}} \right)}} \right\rbrack}}},} & (1.22) \\ {{C_{mn} = {- {\left( {\xi_{m}^{2} - \xi_{n}^{2}} \right)^{- 1}\left\lbrack {F_{mn} - {G_{mn}\xi_{n}^{2}} + {\frac{\partial H}{\partial y}\left( {J_{mn} - {K_{mn}\xi_{n}^{2}}} \right)}} \right\rbrack}}},} & (1.23) \end{matrix}$ and $\begin{matrix} {{B_{nn} = {{- \frac{1}{2}}\left( {E_{nn} + {\frac{\partial H}{\partial x}K_{nn}}} \right)}},} & (1.24) \\ {C_{nn} = {{- \frac{1}{2}}{\left( {G_{nn} + {\frac{\partial H}{\partial y}K_{nn}}} \right).}}} & (1.25) \end{matrix}$

In practice, it is necessary to work with a finite number of coupled mode equations. For this reason, we truncate the coupled mode equations (1.12-1.14) as follows. Fixing a positive integer M, we assume that the expansion coefficients a_(n), b_(n) and c_(n) vanish if n>M. Thus, if n=1, . . . , M then $\begin{matrix} {{{\frac{\partial b_{n}}{\partial x} + \frac{\partial c_{n}}{\partial y} + {\xi_{n}^{2}a_{n}}} = {\sum\limits_{m = 1}^{M}{\left( {{b_{m}B_{mn}} + {c_{m}C_{mn}}} \right)\frac{1}{d_{m}}}}},} & (1.26) \\ {{b_{n} = {\frac{\partial a_{n}}{\partial x} + {\sum\limits_{m = 1}^{M}{a_{m}B_{n\quad m}\frac{1}{d_{m}}}}}},} & (1.27) \\ {c_{n} = {\frac{\partial a_{n}}{\partial y} + {\sum\limits_{m = 1}^{M}{a_{m}C_{n\quad m}{\frac{1}{d_{m}}.}}}}} & (1.28) \end{matrix}$

We can reduce this system of equations to an equivalent system of second-order partial differential equations for just the coefficients a_(n). If n=1, . . . , M then $\begin{matrix} {{{\frac{\partial^{2}a_{n}}{\partial x^{2}} + \frac{\partial^{2}a_{n}}{\partial y^{2}} + {\xi_{n}^{2}a_{n}}} = {\sum\limits_{m = 1}^{M}{\left( {{\alpha_{n\quad m}a_{m}} + {\beta_{n\quad m}\frac{\partial a_{m}}{\partial x}} + {\gamma_{n\quad m}\frac{\partial a_{m}}{\partial y}}} \right)\frac{1}{d_{m}}}}},} & (1.29) \end{matrix}$ where $\begin{matrix} {{a_{n\quad m} = {{\sum\limits_{l = 1}^{M}{\left( {{B_{l\quad m}B_{l\quad n}} + {C_{l\quad m}C_{l\quad n}}} \right)\frac{1}{d_{l}}}} - \frac{\partial B_{n\quad m}}{\partial x} - \frac{\partial C_{n\quad m}}{\partial y}}},} & (1.30) \end{matrix}$

The waveguide considered in the preceding discussion is a slab between two parallel planes. Now let us consider a waveguide that is a shell between two concentric spheres. Let the radius of the outer sphere be R₀, and let the radius of the inner sphere be R₀−L. We pick a spherical coordinate system for this problem anticipating that our analysis of sound propagating in the shell will be restricted to a neighborhood of a great circle E on the outer sphere. We identify E as the equator of this sphere. Let M denote another great circle on the outer sphere that intersects E at right angles and near the source of sound in the waveguide. We identify M as the prime meridian of the outer sphere.

Let P denote a point in the shell. We let τ be the radial distance of P from the common center of the spheres. We let ψ be the angle of latitude of P relative to the equator E. We let φ be the angle of longitude of P relative to the prime meridian M. For definiteness let φ=0 where E intersects M near the source.

It is convenient to introduce a new system of coordinates (x, y, z). x=R₀φ, y=R₀ln(sec ψ+tan ψ), z=r−R₀.  (1.32) Note that y≈R₀ψ if |ψ|<<1. In the new coordinate system the partial differential equation that governs the acoustic pressure p in the waveguide is $\begin{matrix} {{{\frac{\partial\quad}{\partial x}\left( {\frac{1}{\rho}\frac{\partial p}{\partial x}} \right)} + {\frac{\partial\quad}{\partial y}\left( {\frac{1}{\rho}\frac{\partial p}{\partial y}} \right)} + {\cos^{2}\psi\frac{\partial\quad}{\partial z}\left( {\frac{1}{\rho}\frac{r^{2}}{R_{0}^{2}}\frac{\partial p}{\partial x}} \right)} + {\cos^{2}\psi\frac{1}{\rho}\frac{\omega^{2}}{c^{2}}\frac{r^{2}}{R_{0}^{2}}p}} = 0.} & (1.33) \end{matrix}$

We assume that the pressure vanishes on the outer sphere, where z=0, and on the inner sphere, where z=−L: p=0 at z=0 and at z=−L.  (1.34) We also assume that there is a single interface between these spheres where the local sound speed and the local fluid density in the shell can be discontinuous. The radial coordinate of this interface can vary with latitude and with longitude. We describe the interface in the new coordinate system by the equation S(x,y,z)=H(x,y)+z=0,  (1.35) where 0<H(x, y)<L. We can express the continuity conditions at the interface as follows. $\begin{matrix} {{p\quad{and}\quad\frac{1}{\rho}\left( {{\frac{\partial H}{\partial x}\frac{\partial p}{\partial x}} + {\frac{\partial H}{\partial y}\frac{\partial p}{\partial y}} + {\cos^{2}\psi\frac{r^{2}}{R_{0}^{2}}\frac{\partial p}{\partial z}}} \right)}{{{are}\quad{continuous}\quad{across}\quad z} = {- {{H\left( {x,y} \right)}.}}}} & (1.36) \end{matrix}$

The local modes of the spherical shell waveguide are characterized by the eigenfunctions ƒ_(n)(x, y, z) and by the squared wavenumbers ξ_(n) ²(x, y) that satisfy the differential equation $\begin{matrix} {{{{\frac{\partial}{\partial z}\left( {\frac{1}{\rho}\frac{r^{2}}{R_{0}^{2}}\frac{\partial f_{n}}{\partial z}} \right)} + {\frac{1}{\rho}\left( {{\frac{\omega^{2}}{c^{2}}\frac{r^{2}}{R_{0}^{2}}} - \xi_{n}^{2}} \right)f_{n}}} = 0}{{{{if}\quad - L} < z < {0\quad{and}\quad z} \neq {- {H\left( {x,y} \right)}}},}} & (1.37) \end{matrix}$ subject to the conditions that ƒ_(n)=0 at z=−L and at z=0,  (1.38) $\begin{matrix} {{f_{n}\quad{and}\quad\frac{1}{\rho}\frac{\partial f_{n}}{\partial z}\quad{are}\quad{continuous}\quad{across}\quad z} = {- {{H\left( {x,y} \right)}.}}} & (1.39) \end{matrix}$ At each position (x, y) these equations describe a nonsingular selfadjoint characteristic value problem. The squared wavenumbers ξ_(n) ²(x, y) are real and countably infinite in number. We can assume that they are indexed to decrease monotonically to −∞ as the index n increases through the positive integers. We can also assume that σƒ_(n)(x, y, 0)/σz≠0 and that ƒ_(n)(x, y, z) does not vanish identically. Furthermore, we can assume that ƒ_(n)(x, y, z) is real-valued. Since the density ρ is positive, the eigenfunctions are orthogonal in the following sense: $\begin{matrix} {{\int_{- L}^{0}{\frac{1}{\rho}f_{m}f_{n}{\mathbb{d}z}}} = \left\{ \begin{matrix} {{d_{n}\left( {x,y} \right)} > 0} & {{{{if}\quad m} = n},} \\ 0 & {{{if}\quad m} \neq {n.}} \end{matrix} \right.} & (1.40) \end{matrix}$ We assume that the eigenfunctions are normalized so that the signs of the partial derivatives σƒ_(n)(x, y, 0)/σz and the values of the integrals d_(n)(x, y) are independent of the position (x, y).

The functions p, σp/σx and σp/σy have the following eigenfunction expansions: $\begin{matrix} {{p = {{\sum\limits_{n = 1}^{\infty}{a_{n}{f_{n}/d_{n}}\quad{where}\quad a_{n}}} = {\int_{- L}^{0}{\frac{1}{\rho}p\quad f_{n}{\mathbb{d}z}}}}},} & (1.41) \\ {{\frac{\partial p}{\partial x} = {{\sum\limits_{n = 1}^{\infty}{b_{n}f_{n}d_{n}\quad{where}\quad b_{n}}} = {\int_{- L}^{0}{\frac{1}{\rho}\frac{\partial p}{\partial x}f_{n}{\mathbb{d}z}}}}},} & (1.42) \\ {\frac{\partial p}{\partial y} = {{\sum\limits_{n = 1}^{\infty}{c_{n}{f_{n}/d_{n}}\quad{where}\quad c_{n}}} = {\int_{- L}^{0}{\frac{1}{\rho}\frac{\partial p}{\partial y}f_{n}{{\mathbb{d}z}.}}}}} & (1.43) \end{matrix}$ The coefficients a_(n), b_(n) and c_(n) in these expansions are functions of the position (x,y). They are governed by the following system of partial differential equations. $\begin{matrix} {{{\frac{\partial b_{n}}{\partial x} + \frac{\partial c_{n}}{\partial y} + {\xi_{n}^{2}\cos^{2}\psi\quad a_{n}}} = {\sum\limits_{m = 1}^{\infty}{\left( {{b_{m}B_{m\quad n}} + {c_{m}C_{m\quad n}}} \right)\frac{1}{d_{m}}}}},} & (1.44) \\ {{b_{n} = {\frac{\partial a_{n}}{\partial x} + {\sum\limits_{m = 1}^{\infty}{a_{m}B_{m\quad n}\frac{1}{d_{m}}}}}},} & (1.45) \\ {c_{n} = {\frac{\partial a_{n}}{\partial y} + {\sum\limits_{m = 1}^{\infty}{a_{m}C_{n\quad m}{\frac{1}{d_{m}}.}}}}} & (1.46) \end{matrix}$ These equations constitute the fundamental order formulation of coupled mode theory for an irregular three-dimensional waveguide in spherical geometry.

The coupling coefficients B_(mn) and C_(mn) can be computed as follows. Let $\begin{matrix} {g_{n} = {\frac{1}{\rho}\frac{\partial f_{n}}{\partial z}}} & (1.47) \end{matrix}$ and define the following matrix elements: $\begin{matrix} {\quad{{D_{m\quad n} = {\int_{- L}^{0}{\left\lbrack {{\frac{\partial}{\partial x}\left( {\frac{1}{\rho}\frac{w^{2}}{c^{2}}} \right)f_{m}f_{n}} + {\frac{\partial\rho}{\partial x}g_{m}g_{n}}} \right\rbrack\frac{r^{2}}{R_{0}^{2}}{\mathbb{d}z}}}},}} & (1.48) \\ {\quad{{E_{m\quad n} = {\int_{- L}^{0}{\frac{\partial}{\partial x}\left( \frac{1}{\rho} \right)f_{m}f_{n}{\mathbb{d}z}}}},}} & (1.49) \\ {\quad{{F_{m\quad n} = {\int_{- L}^{0}{\left\lbrack {{\frac{\partial}{\partial x}\left( {\frac{1}{\rho}\frac{\omega^{2}}{c^{2}}} \right)f_{m}f_{n}} + {\frac{\partial\rho}{\partial y}g_{m}g_{n}}} \right\rbrack\frac{r^{2}}{R_{0}^{2}}{\mathbb{d}z}}}},}} & (1.50) \\ {\quad{{G_{m\quad n} = {\int_{- L}^{0}{\frac{\partial}{\partial y}\left( \frac{1}{\rho} \right)f_{m}f_{n}{\mathbb{d}z}}}},}} & (1.51) \\ {{J_{m\quad n} = {\left\{ {{\left\lbrack {\frac{1}{\rho}\frac{\omega^{2}}{c^{2}}} \right\rbrack_{{- H} - 0}^{{- H} + 0}{f_{m}\left( {- H} \right)}{f_{n}\left( {- H} \right)}} + {\lbrack\rho\rbrack_{{- H} - 0}^{{- H} + 0}{g_{m}\left( {- H} \right)}{g_{m}\left( {- H} \right)}}} \right\}\quad\left( {1 - {H/R_{0}}} \right)^{2}}},} & (1.52) \\ {\quad{K_{m\quad n} = {\left\lbrack \frac{1}{\rho} \right\rbrack_{{- H} - 0}^{{- H} + 0}{f_{m}\left( {- H} \right)}{{f_{n}\left( {- H} \right)}.}}}} & (1.53) \end{matrix}$ These definitions differ in form from the corresponding definitions for the slab waveguide only by the radial factors in the formulas for the matrix elements D_(mn), F_(mn) and J_(mn).

If m≠n then $\begin{matrix} {{B_{m\quad n} = {- {\left( {\xi_{m}^{2} - \xi_{n}^{2}} \right)^{- 1}\left\lbrack {D_{m\quad n} - {E_{m\quad n}\xi_{n}^{2}} + {\frac{\partial H}{\partial x}\left( {J_{m\quad n} - {K_{m\quad n}\xi_{n}^{2}}} \right)}} \right\rbrack}}},} & (1.54) \\ {{C_{m\quad n} = {- {\left( {\xi_{m}^{2} - \xi_{n}^{2}} \right)^{- 1}\left\lbrack {F_{m\quad n} - {G_{m\quad n}\xi_{n}^{2}} + {\frac{\partial H}{\partial y}\left( {J_{m\quad n} - {K_{m\quad n}\xi_{n}^{2}}} \right)}} \right\rbrack}}},} & (1.55) \end{matrix}$ and $\begin{matrix} {{B_{n\quad n} = {{- \frac{1}{2}}\left( {E_{n\quad n} + {\frac{\partial H}{\partial x}K_{n\quad n}}} \right)}},} & (1.56) \\ {C_{n\quad n} = {{- \frac{1}{2}}{\left( {G_{n\quad n} + {\frac{\partial H}{\partial y}K_{n\quad n}}} \right).}}} & (1.57) \end{matrix}$ These equations are formally identical to the corresponding equations for the slab waveguide.

Let us truncate eqs.(1.44-1.46) by fixing a positive integer M and assuming that the expansion coefficients a_(n), b_(n) and c_(n) vanish if n>M. Thus, if n=1, . . . , M then $\begin{matrix} {{{\frac{\partial b_{n}}{\partial x} + \frac{\partial c_{n}}{\partial y} + {\xi_{n}^{2}\cos^{2}\psi\quad a_{n}}} = {\sum\limits_{m = 1}^{M}{\left( {{b_{m}B_{m\quad n}} + {c_{m}C_{m\quad n}}} \right)\frac{1}{d_{m}}}}},} & (1.58) \\ {{b_{n} = {\frac{\partial a_{n}}{\partial x} + {\sum\limits_{m = 1}^{M}{a_{m}B_{n\quad m}\frac{1}{d_{m}}}}}},} & (1.59) \\ {c_{n} = {\frac{\partial a_{n}}{\partial y} + {\sum\limits_{m = 1}^{M}{a_{m}C_{n\quad m}{\frac{1}{d_{m}}.}}}}} & (1.60) \end{matrix}$ We can reduce this system of equations to an equivalent system of second-order partial differential equations for just the coefficients a_(n). If n=1, . . . , M then $\begin{matrix} {{{\frac{\partial^{2}a_{n}}{\partial x^{2}} + \frac{\partial^{2}a_{n}}{\partial y^{2}} + {\xi_{n}^{2}\cos^{2}\psi\quad a_{n}}} = {\sum\limits_{m = 1}^{M}{\left( {{\alpha_{n\quad m}a_{m}} + {\beta_{n\quad m}\frac{\partial a_{m}}{\partial x}} + {\gamma_{n\quad m}\frac{\partial a_{m}}{\partial y}}} \right)\frac{1}{d_{m}}}}},} & (1.61) \end{matrix}$ where $\begin{matrix} {{\alpha_{n\quad m} = {{\sum\limits_{l = 1}^{M}{\left( {{B_{l\quad m}B_{l\quad n}} + {C_{l\quad m}C_{l\quad n}}} \right)\frac{1}{d_{l}}}} - \frac{\partial B_{n\quad m}}{\partial x} - \frac{\partial C_{n\quad m}}{\partial y}}},} & (1.62) \end{matrix}$ 2. Paraxial approximations for coupled mode equations.

Since the coupled mode equations (1.29,1.61) are elliptic partial differential equations, they can model propagation in all directions at the same time. In many problems of practical interest, however, it is sufficient to consider a single direction of propagation at a time. In slab geometry we can take this direction to be the direction of increasing range, and in spherical geometry we can take it to be the direction of increasing longitude along the equator. In each case, we can formulate parabolic equations from the coupled mode equations to model propagation of sound in the given direction. We shall describe how to formulate these parabolic equations from the coupled mode equations in slab geometry. The approach is identical for the coupled mode equations in spherical geometry.

We let the normalization integral d_(m)=1 for every index m. Hence, the coupled mode equations take the form $\begin{matrix} {{{\frac{\partial^{2}a_{m}}{\partial x^{2}} + \frac{\partial^{2}a_{m}}{\partial y^{2}} + {\xi_{m}^{2}a_{m}}} = {\sum\limits_{n = 1}^{M}\left( {{\alpha_{n\quad m}a_{n}} + {\beta_{m\quad n}\frac{\partial a_{n}}{\partial x}} + {\gamma_{m\quad n}\frac{\partial a_{n}}{\partial y}}} \right)}}{\left( {{m = 1},\ldots\quad,M} \right).}} & (2.1) \end{matrix}$ We introduce a positive reference wavenumber k_(m)(x) for each index m at every range x. We define the corresponding reference phase σ_(m)(x) by $\begin{matrix} {{\vartheta_{m}(x)} = {\int_{0}^{x}{{k_{m}(u)}{{\mathbb{d}u}.}}}} & (2.2) \end{matrix}$ We define the complex amplitude φ_(m)(x, y) of the expansion coefficient a_(m)(x, y) as follows. φ_(m)(x,y)=a_(m)(x,y) exp{−iσ _(m)(x)}.  (2.3) Noting that dσ_(m)(x)/dx=k_(m)(x) and letting k_(m)(x)=dk_(m)(x)/dx, we see that $\begin{matrix} {{\frac{\partial a_{m}}{\partial x} = {\left( {\frac{\partial\varphi_{m}}{\partial x} + {{\mathbb{i}}\quad k_{m}\varphi_{m}}} \right){\mathbb{e}}^{{\mathbb{i}}\quad\vartheta\quad m}}},} & (2.4) \\ {\frac{\partial^{2}a_{m}}{\partial x^{2}} = {\left( {\frac{\partial^{2}\varphi_{m}}{\partial x^{2}} + {2\quad{\mathbb{i}}\quad k_{m}\frac{\partial\varphi_{m}}{\partial x}} + {{\mathbb{i}}{\overset{.}{\quad k}}_{m}\varphi_{m}} - {k_{m}^{2}\varphi_{m}}} \right){{\mathbb{e}}^{{\mathbb{i}}\quad\vartheta\quad m}.}}} & (2.5) \end{matrix}$ Now we assume that the following paraxial approximations are valid: $\begin{matrix} {{\frac{\partial^{2}\varphi_{m}}{\partial x^{2}}} ⪡ {2\quad k_{m}{\frac{\partial\varphi_{m}}{\partial x}}\quad{\left( {{m = 1},\ldots\quad,M} \right).}}} & (2.6) \end{matrix}$ These assumptions allow us to approximate eq.(2.5) by deleting the first term on its right-hand side (RHS). $\begin{matrix} {{\frac{\partial^{2}a_{m}}{\partial x^{2}} = {\left( {{2{\mathbb{i}}\quad k_{m}\frac{\partial\varphi_{m}}{\partial x}} + {{\mathbb{i}}\quad{\overset{.}{k}}_{m}\varphi_{m}} - {k_{m}^{2}\varphi_{m}}} \right){\mathbb{e}}^{{\mathbb{i}}\quad\vartheta\quad m}}},} & (2.7) \end{matrix}$ Combining eqs.(2.3,2.4,2.7) with eq.(2.1), we get the following coupled parabolic equations. $\begin{matrix} {{{2\quad{\mathbb{i}}\quad k_{m}\frac{\partial\varphi_{m}}{\partial x^{2}}} + {{\mathbb{i}}\quad{\overset{.}{k}}_{m}\varphi_{m}} + \frac{\partial^{2}\varphi_{m}}{\partial y^{2}}} = {{\left( {k_{m}^{2} - \xi_{m}^{2}} \right)\varphi_{m}} + {\sum\limits_{n = 1}^{M}{\left\lbrack {{\left( {\alpha_{m\quad n} + {{\mathbb{i}\beta}_{m\quad n}k_{n}}} \right)\varphi_{n}} + {\beta_{m\quad n}\frac{\partial\varphi_{n}}{\partial x}} + {\gamma_{m\quad n}\frac{\partial\varphi_{n}}{\partial y}}} \right\rbrack{\mathbb{e}}^{{\mathbb{i}}{({\vartheta_{m} - \vartheta_{n}})}}}}}} & (2.8) \end{matrix}$ for m=1, . . . , M.

Note the appearance of the partial derivatives σφ_(n)/σx on the RHS of eq.(2.8). Let us rewrite the coupled parabolic equations as follows. $\begin{matrix} {{2\quad{{\mathbb{i}}\left\lbrack {{k_{m}\frac{\partial\varphi_{m}}{\partial x}} + {\frac{i}{2}{\sum\limits_{n = 1}^{M}{\beta_{m\quad n}{\mathbb{e}}^{{- {\mathbb{i}}}\quad{({\vartheta_{m} - \vartheta_{n}})}}\frac{\partial\varphi_{n}}{\partial x}}}}} \right\rbrack}} = {{{- {\mathbb{i}}}\quad k_{m}\varphi_{m}} - \frac{\partial^{2}\varphi_{m}}{\partial y^{2}} + {\left( {k_{m}^{2} - \xi_{m}^{2}} \right)\varphi_{m}} + {\sum\limits_{n = 1}^{M}{\left\lbrack {{\left( {\alpha_{m\quad n} + {{\mathbb{i}}\quad\beta_{m\quad n}k_{n}}} \right)\varphi_{n}} + {\gamma_{m\quad n}\frac{\partial\varphi_{n}}{\partial y}}} \right\rbrack{\mathbb{e}}^{- {{\mathbb{i}}{({\vartheta_{m} - \vartheta_{n}})}}}}}}} & (2.9) \end{matrix}$ for m=1, . . . , M. This system of equations can be solved uniquely for the partial derivatives σφ_(n)/σx if and only if the matrix $M = \quad{(2.10)\left\lbrack \quad\begin{matrix} k_{1} & {\left( {i/2} \right)\beta_{12}e^{- {i{({\vartheta_{1} - \vartheta_{2}})}}}} & \ldots & {\left( {i/2} \right)\beta_{1M}e^{- {i{({\vartheta_{1} - \vartheta_{M}})}}}} \\ {\left( {i/2} \right)\beta_{21}e^{- {i{({\vartheta_{2} - \vartheta_{1}})}}}} & k_{2} & \ldots & {\left( {i/2} \right)\beta_{2\quad M}e^{i{({\vartheta_{2} - \vartheta_{M}})}}} \\ \vdots & \vdots & ⋰ & \vdots \\ {\left( {i/2} \right)\beta_{M1}e^{- {i{({\vartheta_{M} - \vartheta_{1}})}}}} & {\left( {i/2} \right)\beta_{M\quad 2}e^{- {i{({\vartheta_{M} - \vartheta_{2}})}}}} & \ldots & k_{M} \end{matrix} \right\rbrack}$ is nonsingular. (Note that B_(mm)=0 for every m.)

We call M the reference matrix. Since the coupling coefficients B_(mn) are real-valued and skew-symmetric (that is, B_(nm)=−B_(mn) for every m and n), this matrix is Hermitian. Note that M is a function of position (x, y). We shall assume that M is positive definite at every position (x, y). A sufficient condition for M to be positive definite is that this matrix be strictly diagonally dominant. $\begin{matrix} {{{k_{n} > {\left( {1/2} \right){\sum\limits_{m \neq n}{{\beta_{m\quad n}}\quad{for}\quad n}}}} = 1},\ldots\quad,{M.}} & (2.11) \end{matrix}$

There is a useful alternative formulation of the coupled parabolic equations. First, define $\begin{matrix} {\quad{{\psi_{m} = {{k_{m}\varphi_{m}} + {\left( {i/2} \right){\sum\limits_{n = 1}^{M}{\beta_{m\quad n}\varphi_{n}e^{- {i{({\vartheta_{m} - \vartheta_{n}})}}}}}}}},}} & (2.12) \\ {q_{m} = {{\left( {k_{m}^{2} - \xi_{m}^{2}} \right)\varphi_{m}} +}} & (2.13) \\ {\quad{\sum\limits_{n = 1}^{M}\left\lbrack {{\left( {\alpha_{m\quad n} + {\frac{k_{m}}{2\quad k_{m}}\beta_{m\quad n}} + {i\quad k_{m}\beta_{m\quad n}} - \frac{\partial\beta_{m\quad n}}{\partial x}} \right)\varphi_{n}} +} \right.}} & \quad \\ {\left. \quad{{\gamma_{m\quad n}\frac{\partial\varphi_{n}}{\partial y}} + \quad{\frac{i}{2k_{m}}\frac{\partial^{2}}{\partial y^{2}}\left( {\beta_{m\quad n}\varphi_{n}} \right)}} \right\rbrack{e^{{- i}\quad{({\vartheta_{m} - \vartheta_{n}})}}.}} & \quad \end{matrix}$ Now the following equation is equivalent to eq.(2.8). $\begin{matrix} {{{2i\quad k_{m}\frac{\partial\psi_{m}}{\partial x}} - {i\quad k_{m}\psi_{m}} + \frac{\partial^{2}\psi_{m}}{\partial y^{2}}} = {k_{m}{q_{m}.}}} & (2.14) \end{matrix}$ Although eqs.(2.8,2.14) are similar, they differ in one important way: no partial derivative with respect to x of any dependent variable appears on the RHS of eq.(2.14).

Note that we can solve the system of equations (2.12)(m=1, . . . , M) for the complex amplitudes φ_(m) because the reference matrix M is nonsingular. Define the complex column vectors φ=[φ₁, φ₂, . . . , φ_(M)]^(T),  (2.15) ψ=[ψ₁, ψ₂, . . . , ψ_(M)]^(T),  (2.16) where the superscript T indicates transpose. Then, ψ=Mφ, so φ=M⁻¹ψ.

Solutions of the coupled parabolic equations (2.8) obey a conservation law involving two quantities that are closely related to power flow in the waveguide. These quantities are $\begin{matrix} {{J_{P} = {{\sum\limits_{m = 1}^{M}{k_{m}{\varphi_{m}}^{2}}} + {\left( {i/2} \right){\sum\limits_{m,{n = 1}}^{M}{\beta_{m\quad n}{\overset{\_}{\varphi}}_{m}\varphi_{n}e^{- {i{({\vartheta_{m} - \vartheta_{n}})}}}}}}}},} & (2.17) \\ {{K_{P} = {\frac{i}{2}\left\lbrack {{\sum\limits_{m = 1}^{M}\left( {{\varphi_{m}\frac{\partial\varphi_{m}}{\partial y}} - {\frac{\partial\varphi_{m}}{\partial y}\varphi_{m}}} \right)} + {\sum\limits_{m,{n = 1}}^{M}{\gamma_{m\quad n}\varphi_{m}\varphi_{n}e^{- {i{({\vartheta_{m} - \vartheta_{n}})}}}}}} \right\rbrack}},} & (2.18) \end{matrix}$ where {overscore (φ)} is the complex conjugate of φ. The conservation law states that $\begin{matrix} {{\frac{\partial J_{P}}{\partial x} + \frac{\partial K_{P}}{\partial y}} = 0.} & (2.19) \end{matrix}$

Note that J_(p)=φ^(H)Mφ, where φ^(H) 32 |{overscore (φ)}₁, {overscore (φ)}₂, . . . , {overscore (φ)}_(M)| is the Hermitian transpose of φ. If φ≠0, then J_(p)>0 because M is positive definite.

Now let us discuss power flow in the waveguide and show how the quantities J_(p) and K_(p) are related to it. First, consider a vertical strip of infinitesimal width Δy that extends from the bottom of the waveguide at z=−L to the surface of the waveguide at z=0 and that is perpendicular to the x direction. The time-averaged power through the strip is (2ω)⁻¹JΔy, where $\begin{matrix} {J = {\frac{1}{2i}{\int_{- L}^{0}{\frac{1}{\rho}\left( {{\overset{\_}{p}\frac{\partial p}{\partial x}} - {p\frac{\partial\overset{\_}{p}}{\partial x}}} \right){{\mathbb{d}z}.}}}}} & (2.20) \end{matrix}$

Second, consider a vertical strip of infinitesimal width Δx that extends from the bottom of the waveguide to the surface of the waveguide and that is perpendicular to the y direction. The time-averaged power through this strip is (2w)⁻¹KΔx, where $\begin{matrix} {K = {\frac{1}{2i}{\int_{- L}^{0}{\frac{1}{\rho}\left( {{\overset{\_}{p}\frac{\partial p}{\partial y}} - {p\frac{\partial\overset{\_}{p}}{\partial y}}} \right){{\mathbb{d}z}.}}}}} & (2.21) \end{matrix}$

Finally, consider a vertical cylinder that extends from the bottom of the waveguide to the surface of the waveguide and that has a simply-connected cross-section A. Using Green's theorem, we find that the time-averaged power out of this cylinder is $\begin{matrix} {{{\frac{1}{2\quad\omega}{\oint_{\partial A}\left( {{J\quad{\mathbb{d}y}} - {K\quad{\mathbb{d}x}}} \right)}} = {\frac{1}{2\quad\omega}{\int{\int_{A}{\left( {\frac{\partial J}{\partial x} + \frac{\partial K}{\partial y}} \right){\mathbb{d}x}{\mathbb{d}y}}}}}},} & (2.22) \end{matrix}$ where the boundary curve σA has the usual counter-clockwise orientation. Note that the partial derivatives σJ/σx and σK/σy have the physical dimensions of power flux, that is, power/area×time. For this reason, we call J and K power flux integrals.

Suppose that we evaluate these power flux integrals for a solution of the acoustic equation (1.1) that satisfies eq.(1.2) at the boundaries and eq.(1.4) at the interface. We can use these equations and the divergence theorem to show that the integral on the left-hand side (LHS) of eq.(2.22) vanishes. Clearly, the integral on the RHS of this equation vanishes as well. Since the cross-section A is arbitrary, it follows that the integrand of this integral must vanish. Thus, $\begin{matrix} {{\frac{\partial J}{\partial x} + \frac{\partial K}{\partial y}} = 0} & (2.23) \end{matrix}$ for exact solutions of eqs.(1.1,1.2,1.4).

If, however, we evaluate the power flux integrals J and K for an approximate solution of these equations that is constructed from complex amplitude solutions φ_(m) of the coupled parabolic equations (2.8), then we find that $\begin{matrix} {{J = {{\frac{1}{2i}{\sum\limits_{m = 1}^{M}\left( {{{\overset{\_}{\varphi}}_{m}\frac{\partial\varphi_{m}}{\partial x}} - {\varphi_{m}\frac{\partial{\overset{\_}{\varphi}}_{m}}{\partial x}}} \right)}} + J_{P}}},} & (2.24) \\ {K = {K_{P}.}} & (2.25) \end{matrix}$ In this case, eq.(2.19) implies that eq.(2.23) is not valid. Instead, we find that $\begin{matrix} {{\frac{\partial J}{\partial x} + \frac{\partial K}{\partial y}} = {\frac{1}{2i}{\sum\limits_{m = 1}^{M}{\left( {{{\overset{\_}{\varphi}}_{m}\frac{\partial^{2}\varphi_{m}}{\partial x^{2}}} - {\varphi_{m}\frac{\partial^{2}{\overset{\_}{\varphi}}_{m}}{\partial x^{2}}}} \right).}}}} & (2.26) \end{matrix}$ This equation agrees with eq.(2.23) only to the extent that the paraxial approximations (2.6) are satisfied. 3. Adaptive selection of reference wavenumbers to conserve power.

There may be a cross-range W such that we can neglect the acoustic pressure in sections of the waveguide where |y|>W. For example, this can happen if the source produces a bean. In this case, we can assume that a pair of hard vertical walls in the planes through y=±W transversely bound the waveguide. This assumption makes the waveguide a duct, and it imposes the following Neumann boundary conditions on the complex amplitudes φ_(m). $\begin{matrix} {{{\frac{\partial}{\partial y}{\varphi_{m}\left( {x \pm W} \right)}} = {{0\quad{for}\quad m} = 1}},\ldots\quad,{M.}} & (3.1) \end{matrix}$ In this case we also can assume that the coupling coefficients γ_(mn) vanish at the walls. γ_(mn)(x,±W)=0 for m,n=1, . . . , M.  (3.2) Under these assumptions, eqs.(2.18,2.25) imply that the power flux integrals K=K_(p) vanish at the walls. K(x,±W)=K_(p)(x,±W)=0.  (3.3) Thus, if we integrate the conservation law (2.19) from y=−W to y=W, then we find that $\begin{matrix} {{\frac{\partial}{\partial x}{\int_{- W}^{W}{J_{P}{\mathbb{d}y}}}} = 0.} & (3.4) \end{matrix}$ This equation implies that the integral ∫_(−W) ^(W)J_(p)dy is independent of the range.

If eq.(2.23) were valid, then the longitudinal power (2w)⁻¹∫_(−W) ^(W)J dy through the duct would be independent of the range as well. Although eq.(2.23) is not valid for the coupled parabolic equations, we can adaptively select the reference wavenumbers to make the longitudinal power through the duct be independent of the range. If we integrate eq. (2.24) from y=−W to y=W, then we get $\begin{matrix} {{\int_{- W}^{W}{J{\mathbb{d}y}}} = {{\sum\limits_{m = 1}^{M}{\frac{1}{2i}{\int_{- W}^{W}{\left( {{{\overset{\_}{\varphi}}_{m}\frac{\partial\varphi_{m}}{\partial x}} - {\varphi_{m}\frac{\partial{\overset{\_}{\varphi}}_{m}}{\partial x}}} \right){\mathbb{d}y}}}}} + {\int_{- W}^{W}{J_{P}{{\mathbb{d}y}.}}}}} & (3.5) \end{matrix}$ Since the integral ∫_(−W) ^(W)J_(p)dy is independent of the range, this equation implies that the longitudinal power through the duct also will be independent of the range if the sum on the RHS is constant. The best way to accomplish this is to require that each summand vanish, which yields the integral constraints $\begin{matrix} {{{\int_{- W}^{W}{\left( {{{\overset{\_}{\varphi}}_{m}\frac{\partial\varphi_{m}}{\partial x}} - {\varphi_{m}\frac{\partial{\overset{\_}{\varphi}}_{m}}{\partial x}}} \right){\mathbb{d}y}}} = {{0\quad{for}\quad m} = 1}},\ldots\quad,{M.}} & (3.6) \end{matrix}$

These constraints are equivalent to a system of coupled quadratic equations in the reference wavenumbers k_(m). Let k=[k₁, k₂, . . . , k_(M)]^(T). For m=1, . . . , M define the following functions g_(m)(k). $\begin{matrix} \begin{matrix} {{g_{m}(k)} = {{k_{m}^{2}{\int_{- W}^{W}{{\varphi_{m}}^{2}{\mathbb{d}y}}}} + {\sum\limits_{n = 1}^{M}{k_{n}\quad{Re}{\int_{- W}^{W}{\varphi_{m}i\quad\beta_{m\quad n}\varphi_{n}{\mathbb{e}}^{- {{\mathbb{i}}{({\vartheta_{m} - \vartheta_{n}})}}}{\mathbb{d}y}}}}} +}} \\ {{\int_{- W}^{W}{\left( {{\frac{\partial\varphi_{m}}{\partial y}}^{2} - {\xi_{m}^{2}{\varphi_{m}}^{2}}} \right){\mathbb{d}y}}} +} \\ {\sum\limits_{n = 1}^{M}{{Re}\quad{\int_{- W}^{W}{{\varphi_{m}\left( {{\alpha_{m\quad n}\varphi_{n}} + {\beta_{m\quad n}\frac{\partial\varphi_{n}}{\partial x}} + {\gamma_{m\quad n}\frac{\partial\varphi_{n}}{\partial y}}} \right)}{\mathbb{e}}^{- {{\mathbb{i}}{({\vartheta_{m} - \vartheta_{n}})}}}{{\mathbb{d}y}.}}}}} \end{matrix} & (3.7) \end{matrix}$ Then the constraints (3.6) are satisfied if and only if g_(m)(k)=0 for m=1, . . . , M.  (3.8)

It is convenient to write this equation in vector form. Let g(k)=[g₁(k), . . . , g_(M)(k)]^(T). Then eq.(3.8) is equivalent to the equation g(k)=0. In practice we can use Newton's method to solve this equation iteratively. The Jacobian matrix σg/σk that is associated with g(k) has the elements $\begin{matrix} {\frac{\partial g_{m}}{\partial k_{n}} = \left\{ \begin{matrix} {2k_{m}{\int_{- W}^{W}{{\varphi_{m}}^{2}{\mathbb{d}y}}}} & {{{{if}\quad m} = n},\left( {{{since}\quad\beta_{m\quad m}} = 0} \right)} \\ {{Re}{\int_{- W}^{W}{{\overset{\_}{\varphi}}_{m}i\quad\beta_{mn}\varphi_{n}{\mathbb{e}}^{- {{\mathbb{i}}{({\theta_{m} - \theta_{n}})}}}{\mathbb{d}y}}}} & {{{if}\quad m} \neq {n.}} \end{matrix} \right.} & (3.9) \end{matrix}$ Given an estimate k₀ of the solution vector k, we can compute a correction Δk to it as the solution of the linear equation $\begin{matrix} {\left. \frac{\partial g}{\partial k} \middle| {}_{k_{0}}{{\cdot \Delta}\quad k} \right. = {- {{g\left( k_{0} \right)}.}}} & (3.10) \end{matrix}$ Then we replace k₀ with k₀+Δk. This completes one iteration of Newton's method. The Jacobian matrix σg/σk is real symmetric. If no amplitude φ_(m) vanishes identically and if the reference matrix M is positive definite at every position, then the Jacobian matrix σg/σk also is positive definite.

The solutions of eq.(3.8) are optimal, but they are expensive to compute. There are suboptimal ways to estimate these reference wavenumbers cheaply that can be useful. For example, if we neglect the coupling coefficients in eq.(3.7), then we get the following approximate solution of eq.(3.8). $\begin{matrix} {{k_{m}^{2} = {\int_{- W}^{W}{\left( {{\xi_{m}^{2}{\varphi_{m}}^{2}} - {\frac{\partial\varphi_{m}}{\partial y}}^{2}} \right){{\mathbb{d}y}/{\int_{- W}^{W}{{\varphi_{m}}^{2}{\mathbb{d}y}}}}}}}\quad{{{{for}\quad m} = 1},\ldots\quad,{M.}}} & (3.11) \end{matrix}$ 4. Exact integrals of the coupled parabolic equations.

We can use Fourier cosine expansions to integrate the alternative coupled parabolic equations (2.14) subject to Neumann boundary conditions at the boundaries y=±W, $\begin{matrix} {{{\frac{\partial}{\partial y}{\varphi_{m}\left( {x,{\pm W}} \right)}} = {{0\quad{for}\quad m} = 1}},\ldots\quad,{M.}} & (4.1) \end{matrix}$ It is convenient to introduce the following operator notation. Let g(y) be a continuously differentiable function over the interval −W≦y≦W such that dg(±W)/dy=0. The Fourier cosine coefficients ĝ of g are defined by the integrals $\begin{matrix} {{{\hat{g}}_{n} = {{\hat{\varepsilon}}_{n}{\int_{0}^{2W}{{g\left( {s - W} \right)}\cos\left\{ {\left( {n\quad{\pi/2}W} \right)s} \right\}{{\mathbb{d}s}/W}}}}}\quad{{{{for}\quad n} = 0},1,2,\ldots\quad,}} & (4.2) \end{matrix}$ where {circumflex over (ε)}₀=½ and {circumflex over (ε)}_(n)=1 for n=1, 2, . . . . Then the Fourier cosine expansion $\begin{matrix} {{g\left( {s - W} \right)} = {\sum\limits_{n = 0}^{\infty}{{\hat{g}}_{n}\cos\left\{ {\left( {n\quad{\pi/2}W} \right)s} \right\}\quad\left( {0 \leq s \leq {2W}} \right)}}} & (4.3) \end{matrix}$ is uniformly convergent. Let ĝ denote the sequence {ĝ₀, ĝ₁, . . . , ĝ_(n) . . . }. The function g and the sequence ĝ uniquely determine each other. We define the operators ℑ and ℑ⁻¹ by the relations ℑg=ĝ, ℑ⁻¹ĝ=g.  (4.4)

Next, we introduce the spectral parameter η. This parameter takes on the discrete values η=nπ/2W for n=0,1,2, . . .  (4.5) Given a function h(η) and a sequence ℑg, we define the new sequence h(η)ℑg as follows. h(η)ℑg={h(0)ĝ₀, h(π/2W)ĝ₁, . . . , h(ηπ/2W)ĝ_(n), . . . }.  (4.6)

The integral form of eq.(2.14) is ψ m ⁡ ( x j + 1 ) = ( k m ⁡ ( x j + 1 ) k m ⁡ ( x j ) ) 1 2 ⁢ - 1 ⁢ ( exp ⁢ { - ( i ⁢   ⁢ η 2 / 2 ) ⁢ ∫ x j x j + 1 ⁢ ⅆ u k m ⁡ ( u ) } ⁢ ⁢   ⁢ ψ m ⁡ ( x j ) ) - ( i / 2 ) ⁢ ∫ x j x j + 1 ⁢ ( k m ⁡ ( x j + 1 ) k m ⁡ ( x ) ) 1 2 ⁢ - 1 ⁢   ( exp ⁢ { - ( i ⁢   ⁢ η 2 / 2 ) ⁢ ∫ x j x j + 1 ⁢ ⅆ u k m ⁡ ( u ) } ⁢ ⁢   ⁢ q m ⁡ ( x ) ) ⁢ ⅆ x . ( 4.7 ) This equation, which holds for m=1, . . . , M, is the basis of our method for integrating the coupled parabolic equations numerically. 5. Discrete approximation of the exact integrals.

We shall introduce three numerical approximations in eq.(4.7) that lead to a system of linear equations for the values of the complex amplitudes (m at discrete mesh points in the x, y plane. First, we use the trapezoidal rule to approximate the integral from x=x_(j) to x=x_(j+1) on the RHS of this equation. Letting Δx=x_(j+1)−x_(j) and neglecting terms of ο((Δx)³), we get ψ m ⁡ ( x j + 1 ) + ( i ⁢   ⁢ Δ ⁢   ⁢ x / 4 ) ⁢ q m ⁡ ( x j + 1 ) = ( k m ⁡ ( x j + 1 ) k m ⁡ ( x j ) ) 1 2 ×   - 1 ⁢ ( exp ⁢ { - ( i ⁢   ⁢ η 2 / 2 ) ⁢ ∫ x j x j + 1 ⁢ ⅆ u k m ⁡ ( u ) } ⁢ ⁡ [ ψ m ⁡ ( x j ) + ( i ⁢   ⁢ Δ ⁢   ⁢ x / 4 ) ⁢ q m ⁡ ( x j ) ] ) . ( 5.1 )

We also can use the trapezoidal rule to approximate the argument of the exponential function in eq.(5.1). Neglecting terms of ο((Δx)³), we get $\begin{matrix} {{\int_{x_{j}}^{x_{j + 1}}\frac{\mathbb{d}u}{k_{m}(u)}} = {\left( {{k_{m}^{- 1}\left( x_{j} \right)} + {k_{m}^{- 1}\left( x_{j + 1} \right)}} \right)\Delta\quad{x/2.}}} & (5.2) \end{matrix}$ We can approximate the phases σ_(m)(x_(j+1)) the same way. Neglecting terms of ο((Δx)³), we get $\begin{matrix} \begin{matrix} {{\vartheta_{m}\left( x_{j + 1} \right)} = {{\vartheta_{m}\left( x_{j} \right)} + {\int_{x_{j}}^{x_{j + 1}}{{k_{m}(u)}{\mathbb{d}u}}}}} \\ {= {{\vartheta_{m}\left( x_{j} \right)} + {\left( {{k_{m}\left( x_{j} \right)} + {k_{m}\left( x_{j + 1} \right)}} \right)\Delta\quad{x/2.}}}} \end{matrix} & (5.3) \end{matrix}$

The second numerical approximation that we introduce is to discretize the operators ℑ and ℑ⁻¹ using discrete cosine transforms. Fix a positive integer N, and let g={g₀, g₁, . . . , g_(N)) be a sequence of N+1 complex numbers. Define the discrete Fourier cosine coefficients {overscore (g)}_(n) as follows. $\begin{matrix} {{{\overset{\_}{g}}_{n} = {{{{\overset{\_}{\varepsilon}}_{n}\left( {2/N} \right)}{\sum\limits_{k = 0}^{N}{{\overset{\_}{\varepsilon}}_{k}g_{k}{\cos\left( {\pi\quad{{kn}/N}} \right)}\quad{for}\quad n}}} = 0}},1,\ldots\quad,N,{{{where}\quad{\overset{\_}{\varepsilon}}_{0}} = {{\overset{\_}{\varepsilon}}_{N} = {{\frac{1}{2}\quad{and}\quad{\overset{\_}{\varepsilon}}_{n}} = {{1\quad{for}\quad n} = 1}}}},\ldots\quad,{N - {1.\quad{Now}}}} & (5.4) \end{matrix}$  for n=0, 1, . . . , N,  (5.4) where {overscore (ε)}₀={overscore (ε)}_(N)=½ and {overscore (ε)}_(n)=1 for n=1, . . . , N−1. Now $\begin{matrix} {{g_{k} = {{\sum\limits_{n = 0}^{N}{{\overset{\_}{g}}_{n}{\cos\left( {\pi\quad{{kn}/N}} \right)}\quad{for}\quad k}} = 0}},1,\ldots\quad,N} & (5.5) \end{matrix}$  for k=0, 1, . . . , N  (5.5) if and only if the coefficients {overscore (g)}_(n) are given by eq.(5.4). Let {overscore (g)}={{overscore (g)}₀, {overscore (g)}₁ , . . . , {overscore (g)}_(N)} denote the sequence of these coefficients. We define the operators ℑ_(N) and ℑ_(N) ⁻¹ by the relations ℑ_(N)g={overscore (g)}, ℑ_(N) ⁻¹{overscore (g)}=g  (5.6)

We replace the operators ℑ and ℑ⁻¹ in eq.(5.1) by the operators ℑ_(N) and ℑ_(n) ⁻¹. Since these operators transform sequences to sequences, however, we must also replace the functions ψ_(m)±(iΔx/4)q_(m) in this equation with the sequences of their values at the points y_(k), which we write as {ψ_(m)(y_(k))±(iΔx/4)q_(m)(y_(k))}_(k=0) ^(N). Under these approximations, eq.(5.1) becomes { ψ m ⁡ ( x j + 1 , y k ) + ( i ⁢   ⁢ Δ ⁢   ⁢x / 4 ) ⁢ q m ⁡ ( x j + 1 , y k ) } k = 0 N = ( k m ⁡ ( x j + 1 ) k m ⁡ ( x j ) ) 1 2 ×   N - 1 ⁢ ( exp ⁢ { - ( i ⁢   ⁢ η 2 / 2 ) ⁢ ∫ x j x j + 1 ⁢ ⅆ u k m ⁡ ( u ) } ⁢ N ⁢ { ψ m ⁡ ( x j , y l ) - ( i ⁢   ⁢ Δ ⁢   ⁢ x / 4 ) ⁢ q m ⁡ ( x j , y l ) } l = 0 N ) . ( 5.7 ) Note that the spectral parameter η in this equation is bounded because η=ηπ/2W for n=0, 1, . . . , N. Hence, 0≦η≦π/Δy.

We must approximate the partial derivatives σφ_(n)/σy and σ²(β_(mn)φ_(n))/σy² that appear in eq.(2.13) for q_(m) before we can evaluate the sequences in eq.(5.7). Since the variable y in this equation is restricted to the uniformly spaced points y_(k), our third approximation is to replace the partial derivatives at these points with central differences. At the lowest order of approximation, we use the three-point formulas $\begin{matrix} {{{\frac{\partial\quad}{\partial y}{\varphi_{n}\left( y_{k} \right)}} = {{\frac{1}{2\Delta\quad y}\left\lbrack {{\varphi_{n}\left( y_{k + 1} \right)} - {\varphi_{n}\left( y_{k - 1} \right)}} \right\rbrack} + {O\left( \left( {\Delta\quad y} \right)^{2} \right)}}},} & (5.8) \\ {{\frac{\partial^{2}\quad}{\partial y^{2}}\left( {\beta_{mn}\varphi_{n}} \right)\left( y_{k} \right)} = {{\frac{1}{\left( {\Delta\quad y} \right)^{2}}\left\lbrack {{\left( {\beta_{mn}\varphi_{n}} \right)\left( y_{k + 1} \right)} - {2\left( {\beta_{mn}\varphi_{n}} \right)\left( y_{k} \right)} + {\left( {\beta_{mn}\varphi_{n}} \right)\left( y_{k - 1} \right)}} \right\rbrack} + {{O\left( \left( {\Delta\quad y} \right)^{2} \right)}.}}} & (5.9) \end{matrix}$ Now we combine eqs.(5.8,5.9) with eq.(2.13). Rearranging terms in the result, we find that $\begin{matrix} {{q_{m}\left( y_{k} \right)} = {{\left( {k_{m}^{2} - {\xi_{m}^{2}\left( y_{k} \right)}} \right){\varphi_{m}\left( y_{k} \right)}} + {\sum\limits_{n = 1}^{M}{\left\lbrack {{\alpha_{mn}\left( y_{k} \right)} + {\left( {\frac{{\overset{.}{k}}_{m}}{2k_{m}} + {ik}_{m} - \frac{i}{{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right){\beta_{mn}\left( y_{k} \right)}} - {\frac{\partial\quad}{\partial x}{\beta_{mn}\left( y_{k} \right)}}} \right\rbrack e^{- {i{({\theta_{m} - \theta_{n}})}}}{\varphi_{n}\left( y_{k} \right)}}} + {\sum\limits_{n = 1}^{M}{\left( {{\frac{1}{2\Delta\quad y}{\gamma_{mn}\left( y_{k} \right)}} + {\frac{i}{2{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( y_{k + 1} \right)}}} \right)e^{- {i{({\theta_{m} - \theta_{n}})}}}{\varphi_{n}\left( y_{k + 1} \right)}}} + {\sum\limits_{n = 1}^{M}{\left( {{{- \frac{1}{2\Delta\quad y}}{\gamma_{mn}\left( y_{k} \right)}} + {\frac{i}{2{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( y_{k - 1} \right)}}} \right)e^{- {i{({\theta_{m} - \theta_{n}})}}}{\varphi_{n}\left( y_{k - 1} \right)}}} + {{O\left( \left( {\Delta\quad y} \right)^{2} \right)}.}}} & (5.10) \end{matrix}$

Therefore, neglecting terms of ο(Δx(Δy)²), we find that $\begin{matrix} {{{\psi_{m}\left( {x_{j + 1},y_{k}} \right)} + {\left( {i\quad\Delta\quad{x/4}} \right){q_{m}\left( {x_{j + 1},y_{k}} \right)}}} = {\quad{{{\left\lbrack {k_{m} + {\left( {i\quad\Delta\quad{x/4}} \right)\left( {k_{m}^{2} - {\xi_{m}^{2}\left( {x_{j + 1},y_{k}} \right)}} \right)}} \right\rbrack{\varphi_{m}\left( {x_{j + 1},y_{k}} \right)}} + {\left( {\frac{i}{2} + \frac{i{\overset{.}{k}}_{m}\Delta\quad x}{8k_{m}} - \frac{k_{m}\Delta\quad x}{4} + \frac{\Delta\quad x}{4{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right){\sum\limits_{n = 1}^{M}{{\beta_{mn}\left( {x_{j + 1},y_{k}} \right)}e^{- {i{({\theta_{m} - \theta_{n}})}}}{\varphi_{n}\left( {x_{j + 1},y_{k}} \right)}}}} + {\left( {i\quad\Delta\quad{x/4}} \right){\sum\limits_{n = 1}^{M}{\left( {{\alpha_{mn}\left( {x_{j + 1},y_{k}} \right)} - {\frac{\partial\quad}{\partial x}{\beta_{mn}\left( {x_{j + 1},y_{k}} \right)}}} \right)e^{- {i{({\theta_{m} - \theta_{n}})}}}{\varphi_{n}\left( {x_{j + 1},y_{k}} \right)}}}} + {\sum\limits_{n = 1}^{M}{\left( {{\frac{i\quad\Delta\quad x}{8\Delta\quad y}{\gamma_{mn}\left( {x_{j + 1},y_{k}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j + 1},y_{k + 1}} \right)}}} \right)e^{- {i{({\theta_{m} - \theta_{n}})}}}{\varphi_{n}\left( {x_{j + 1},y_{k + 1}} \right)}}} + {\sum\limits_{n = 1}^{M}{\left( {{{- \frac{i\quad\Delta\quad x}{8\Delta\quad y}}{\gamma_{mn}\left( {x_{j + 1},y_{k}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j + 1},y_{k - 1}} \right)}}} \right)e^{- {i{({\theta_{m} - \theta_{n}})}}}{\varphi_{n}\left( {x_{j + 1},y_{k - 1}} \right)}}}},}}} & (5.11) \end{matrix}$ where k_(n)=k_(n)(x_(j+1)) and σ_(n)=σ_(n)(x_(j+1)) for n=1, . . . , M.

Similarly, neglecting terms of ο(Δx(Δy)²), we find that $\begin{matrix} {{{{\psi_{m}\left( {x_{j},y_{l}} \right)} - {\left( {i\quad\Delta\quad{x/4}} \right){q_{m}\left( {x_{j},y_{l}} \right)}}} = {{\left\lbrack \quad{k_{m} - {\left( {i\quad\Delta\quad{x/4}} \right)\left( {k_{m}^{2} - {\xi_{m}^{2}\left( {x_{j},y_{l}} \right)}} \right)}} \right\rbrack\quad{\varphi_{m}\left( {x_{j},y_{l}} \right)}} + {\left( {\frac{i}{2} - \frac{{ik}_{m}\Delta\quad x}{8k_{m}} + \frac{k_{m}\Delta\quad x}{4} - \frac{\Delta\quad x}{4{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right)\quad{\sum\limits_{n = 1}^{M}{{\beta_{mn}\left( {x_{j},y_{l}} \right)}\quad e^{- {i{({\theta_{m} - \theta_{n}})}}}\quad{\varphi_{n}\left( {x_{j},y_{l}} \right)}}}} - {\left( {i\quad\Delta\quad{x/4}} \right){\sum\limits_{n = 1}^{M}{\left( {{\alpha_{mn}\left( {x_{j},y_{l}} \right)} - {\frac{\partial\quad}{\partial x}{\beta_{mn}\left( {x_{j},y_{l}} \right)}}} \right)e^{- {i{({\theta_{m} - \theta_{n}})}}}{\varphi_{n}\left( {x_{j},y_{l}} \right)}}}} + {\sum\limits_{n = 1}^{M}{\left( {{{- \frac{i\quad\Delta\quad x}{8\Delta\quad y}}\gamma_{mn}\left( {x_{j},y_{l}} \right)} + {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j},y_{l + 1}} \right)}}} \right)\quad e^{- {i{({\theta_{m} - \theta_{n}})}}}\quad{\varphi_{n}\left( \quad{x_{j},y_{l + 1}} \right)}}} + {\sum\limits_{n = 1}^{M}{\left( {{{- \frac{i\quad\Delta\quad x}{8\Delta\quad y}}\gamma_{mn}\left( {x_{j},y_{l}} \right)} + {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}\quad{\beta_{mn}\left( \quad{x_{j},\quad y_{l - 1}} \right)}}} \right)\quad e^{- {i{({\theta_{m} - \theta_{n}})}}}\quad{\varphi_{n}\left( \quad{x_{j},y_{l - 1}} \right)}}}}},} & (5.12) \end{matrix}$ where k_(n)=k_(n)(x_(j)) and σ_(n)=σ_(n)(x_(j)) for n=1, . . . , M.

We need additional information to interpret these equations at the boundaries y=±W. Let us require that eq.(4.1) be consistent with eq (3.1). It follows that for k=0, N and for all x_(j), φ_(n)(x_(j),y_(k−1))+φ_(n)(x_(j),y_(k+1)) for n=1, . . . ,M,  (5.13) β_(mn)(x_(j),y_(k−1))=β_(mn)(x_(j),y_(k+1)) for m,n=1, . . . ,M.  (5.14)

Equations (5.7,5.11), subject to eqs (5.13,5.14), determine a system of linear equations for the values φ_(m)(x_(j+1), y_(k)) at a range step x_(j+1), which in turn depend linearly on the values φ_(m)(x_(j), y₁) at the previous range step x_(j) through eq.(5.12). In order to solve these equations numerically, it is helpful to put them in matrix form. First, we write eq.(5.7) as follows. { Ψ m ⁡ ( x j + 1 - 0 , y k ) } k = 0 N = ( k m ⁡ ( x j + 1 ) k m ⁡ ( x j ) ) 1 2 × N - 1 ⁢   ⁢ ( exp ⁢ { - ( i ⁢   ⁢ η 2 / 2 ) ⁢ ∫ x i x i + 1 ⁢ ⅆ u k m ⁡ ( u ) } ⁢ N ⁢ { Ψ m ⁡ ( x j + 0 , y l ) } l = 0 N ) , ( 5.15 ) where Ψ_(m)(x_(j+1)−0,y_(k))=ψ_(m)(x_(j+1),y_(k))+(iΔx/4)q_(m)(x_(j+1),y_(k)),  (5.16) Ψ_(m)(x_(j)+0,y₁)=ψ_(m)(x_(j),y₁)−(iΔx/4)q_(m)(x_(j),y₁).  (5.17) Second, for k=0,1, . . . , N we define the column vectors φ_(k)=[φ₁(x_(j+1),y_(k)), . . . , φ_(M)(x_(j+1),y_(k))]^(T),  (5.18) Ψ_(k) ⁻=[Ψ₁(x_(j+1)−ο,y_(k)), . . . , Ψ_(M)(x_(j+1)−ο, y_(k))]^(T).  (5.19) Third, for k=0,1, . . . , N we define the M×M matrices A_(kk), A_(kk+1) and A_(kk−1) as follows. If k_(n)=k_(n)(x_(j+1)) and σ_(n)=σ_(n)(x_(j+1)) for n=1, . . . , M, then $\begin{matrix} {\left( A_{kk} \right)_{mn} = \left\{ \begin{matrix} {{k_{n} + {\left( {i\quad\Delta\quad{x/4}} \right)\left( {k_{n}^{2} - {\xi_{n}^{2}\left( {x_{j + 1},y_{k}} \right)} + {\alpha_{nn}\left( {x_{j + 1},y_{k}} \right)}} \right)}},} & \left( {{{if}\quad m} = n} \right) \\ \left\lbrack {\left( {\frac{i}{2} + \frac{i{\overset{.}{k}}_{m}\Delta\quad x}{8k_{m}} - \frac{k_{m}\Delta\quad x}{4} + \frac{\Delta\quad x}{4{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right){\beta_{mn}\left( {x_{j + 1},y_{k}} \right)}} \right. & \quad \\ {{{+ \left( {i\quad\Delta\quad{x/4}} \right)}\left( {{\alpha_{mn}\left( {x_{j + 1},y_{k}} \right)} - {\frac{\partial\quad}{\partial x}{\beta_{mn}\left( {x_{j + 1},y_{k}} \right)}}} \right)e^{- {i{({\theta_{m} - \theta_{n}})}}}},} & \left( {{{if}\quad m} \neq n} \right) \end{matrix} \right.} & (5.20) \\ {\left( A_{{kk} \pm 1} \right)_{mn} = {\left( {{{\pm \frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}}{\gamma_{mn}\left( {x_{j + 1},y_{k}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j + 1},y_{k \pm 1}} \right)}}} \right){e^{- {i{({\theta_{m} - \theta_{n}})}}}.}}} & (5.21) \end{matrix}$ Finally, we define the M×M matrices B₁ and B_(N−1) as follows. For k=1, N−1, $\begin{matrix} {\left( B_{k} \right)_{mn} = {{- \frac{\Delta\quad x}{4{k_{m}\left( {\Delta\quad y} \right)}^{2}}}{\beta_{mn}\left( {x_{j + 1},y_{k}} \right)}{e^{- {i{({\theta_{m} - \theta_{n}})}}}.}}} & (5.22) \end{matrix}$ Note that B₁=A_(ο−1)+A_(ο1) and B_(N−1)=A_(NN−1)+A_(NN+1). Now eqs.(5.11,5.16), subject to eqs.(5.13,5.14), are equivalent to the following matrix equation. $\begin{matrix} {{\begin{bmatrix} A_{00} & B_{1} & \quad & \quad & \quad \\ A_{10} & A_{11} & A_{12} & \quad & \quad \\ \quad & ⋰ & ⋰ & ⋰ & \quad \\ \quad & \quad & A_{N - {1N} - 2} & A_{N - {1N} - 1} & A_{N - {1N}} \\ \quad & \quad & \quad & B_{N - 1} & A_{NN} \end{bmatrix} \cdot \begin{bmatrix} \varphi_{0} \\ \varphi_{1} \\ \vdots \\ \varphi_{N - 1} \\ \varphi_{N} \end{bmatrix}} = {\begin{bmatrix} \Psi_{0}^{-} \\ \Psi_{1}^{-} \\ \vdots \\ \Psi_{N - 1}^{-} \\ \Psi_{N}^{-} \end{bmatrix}.}} & (5.23) \end{matrix}$

We recommend solving this equation numerically by Gaussian elimination with partial pivoting (GEPP) followed by one round of iterative refinement. The block tri-diagonal system matrix on the LHS of eq.(5.23) is also a band matrix of order M(N+1) that has 2M−1 lower diagonals and 2M−1 upper diagonals. Since M<<N in practice, it is most efficient to do the GEPP computations with subroutines that take advantage of band, storage.

Let us write eq.(5.23) as A·Φ=Ψ⁻, and let Φ₀ be the numerical solution of this equation computed in finite precision arithmetic by GEPP. Iterative refinement is a technique to improve the accuracy of this solution. One round of this technique has three steps.

(1) Compute the residual τ₀=Ψ⁻−A·Φ₀.

(2) Using the LU decomposition of A from the original solution of eq. (5.23) by GEPP, solve the equation A·δ₀=τ₀ for the correction δ₀.

(3) Refine the solution: Φ_(t)=Φ₀+δ₀.

The computations in these steps are done at the same arithmetic precision as the original solution.

The system matrix on the LHS of eq.(5.23) is nonsingular if it is strictly diagonally dominant by columns. We see from eqs.(5.20,5.21) that this condition requires the following strict inequalities to be satisfied for the block indices k=2, . . . , N−2: $\begin{matrix} {{{{{k_{n} + {\left( {i\quad\Delta\quad{x/4}} \right)\left( {k_{n}^{2} - {\xi_{n}^{2}\left( {x_{j + 1},y_{k}} \right)} + {a_{nn}\left( {x_{j + 1},y_{k}} \right)}} \right)}}} > {{\sum\limits_{\underset{({m \neq n})}{m = 1}}^{M}{{{\left( {\frac{i}{2} + \frac{{ik}_{m}\Delta\quad x}{8k_{m}} - \frac{k_{m}\Delta\quad x}{4} + \frac{\Delta\quad x}{4{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right){\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}} + {\left( {i\quad\Delta\quad{x/4}} \right)\left( {{a_{m\quad n}\left( {x_{j + 1},y_{k}} \right)} - {\frac{\partial\quad}{\partial x}{\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}}} \right)}}}} + {\sum\limits_{m = 1}^{M}{{{\frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}{\gamma_{m\quad n}\left( {x_{j + 1},y_{k - 1}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}}}}} + {\sum\limits_{m = 1}^{M}{{{{\frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}{\gamma_{m\quad n}\left( {x_{j + 1},y_{k + 1}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}}}}\quad{for}\quad n}}}} = 1},\ldots\quad,{M.}} & (5.24) \end{matrix}$ These inequalities are modified slightly for the remaining indices k=0, 1, N−1, N. This condition is more restrictive than the similar condition (2.11) for the reference matrix M to be positive definite at the points y_(k). If this condition is satisfied, then Gaussian elimination requires no pivoting and is stable.

Note the factors k_(m)Δx, Δx/Δy and Δx/k_(m)(Δy)² that appear in this equation. We recommend choosing the increments Δx and Δy so that k_(m)Δx≦ο(1) and k_(m)Δy≧ο(1). These restrictions imply that the ratios Δx/Δy and Δr/k_(m)(Δy)² are at most ο(1). These restrictions also limit the errors that we make by applying the trapezoidal rule to eq.(4.7) to obtain eq.(5.1).

Finally, we note that higher order central difference formulas can be used in place of eqs.(5.8,5.9). For example, at the next order of approximation, we would use five-point formulas. They lead to a matrix formulation like eq.(5.23) in which the system matrix is block penta-diagonal.

6. Discrete transparent boundary conditions.

The most simple standard parabolic equation is $\begin{matrix} {{{{2{ik}_{0}\frac{\partial\varphi}{\partial x}} + \frac{\partial^{2}\varphi}{\partial y^{2}}} = 0},} & (6.1) \end{matrix}$ where the wavenumber k₀ is a positive constant. Let us assume

(1) the solution φ(x, y) is defined in a region of the x, y plane that includes the upper right quadrant where x≧0 and y≧0,

(2) the solution φ(x, y) and its partial derivatives are bounded in this quadrant,

(3) φ(0, y)=0 for y≧0. We can translate these assumptions into a transparent boundary condition that the solution φ(z, y) satisfies on the right-half x-axis (x>0, y=0). This condition makes the half-axis transparent in the sense that the wavefield described by the solution φ(x, y) propagates across this half-axis into the upper right quadrant without reflection. $\begin{matrix} {{\varphi\left( {x,0} \right)} = {{- {e^{i\quad{\pi/4}}\left( {2\pi\quad k_{0}} \right)}^{{- 1}/2}}{\int_{0}^{x}{\frac{\partial\quad}{\partial y}{\varphi\left( {u,0} \right)}\left( {x - u} \right)^{{- 1}/2}{{\mathbb{d}u}.}}}}} & (6.2) \end{matrix}$ If we change the sign of the partial derivative σ/σy in this equation, then we get a new condition that makes the right-half x-axis transparent in the opposite sense. Now the wavefield described by the solution φ(x, y) propagates across this half-axis into the lower right quadrant without reflection. $\begin{matrix} {{\varphi\left( {x,0} \right)} = {{e^{i\quad{\pi/4}}\left( {2\pi\quad k_{0}} \right)}^{{- 1}/2}{\int_{0}^{x}{\frac{\partial\quad}{\partial y}{\varphi\left( {u,0} \right)}\left( {x - u} \right)^{{- 1}/2}{{\mathbb{d}u}.}}}}} & (6.3) \end{matrix}$ Let us replace the reference wavenumber k₀ in eq.(6.1) with a function k(x)>0, $\begin{matrix} {{{2{{ik}(x)}\frac{\partial\varphi}{\partial x}} + \frac{\partial^{2}\varphi}{\partial y^{2}}} = 0.} & (6.4) \end{matrix}$ The change of variable $\begin{matrix} {{\tau(x)} = {\int_{0}^{x}\frac{\mathbb{d}u}{k(u)}}} & (6.5) \end{matrix}$ transforms this equation into a standard parabolic equation, $\begin{matrix} {{{2i\frac{\partial\varphi}{\partial\tau}} + \frac{\partial^{2}\varphi}{\partial y^{2}}} = 0.} & (6.6) \end{matrix}$ Let {x₀,x₁, . . . ,xj, . . . } be an increasing sequence of discrete range steps, with x₀=0. Define ψ_(j)(y)=φ(x, y) and τ_(j)=τ(x_(j)) for j=0,1, . . . . Using the trapezoidal rule to integrate eq.(6.6) from τ_(j) to τ_(j+1) and neglecting the ο((τ_(j+1)−τ_(j))³) remainder, we find that $\begin{matrix} {{{2{i\left\lbrack {{\psi_{j + 1}(y)} - {\psi_{j}(y)}} \right\rbrack}} + {\frac{1}{2}{\left( {\tau_{j + 1} - \tau_{j}} \right)\left\lbrack {{\frac{\partial^{2}}{\partial y^{2}}{\psi_{j}(y)}} + {\frac{\partial^{2}}{\partial y^{2}}{\psi_{j + 1}(y)}}} \right\rbrack}}} = 0.} & (6.7) \end{matrix}$ Similarly, $\begin{matrix} {{\tau_{j + 1} - \tau_{j}} = {{\int_{x_{j}}^{x_{j + 1}}\frac{d\quad u}{k(u)}} = {\frac{1}{2}{{\left( {x_{j + 1} - x_{j}} \right)\left\lbrack {{k^{- 1}\left( x_{j} \right)} + {k^{- 1}\left( x_{j + 1} \right)}} \right\rbrack}.}}}} & (6.8) \end{matrix}$ where we evaluated the integral by the trapezoidal rule and dropped the ο((x_(j+1)−x_(j))³) remainder. We get a semi-discrete parabolic equation from eqs.(6.7,6.8). $\begin{matrix} {{{\frac{\partial^{2}}{\partial y^{2}}{\psi_{j}(y)}} + {\frac{\partial^{2}}{\partial y^{2}}{\psi_{j + 1}(y)}} + {i{\frac{8}{\left( {x_{j + 1} - x_{j}} \right)\left\lbrack {{k^{- 1}\left( x_{j} \right)} + {k^{- 1}\left( x_{j + 1} \right)}} \right\rbrack}\left\lbrack {{\psi_{j + 1}(y)} - {\psi_{j}(y)}} \right\rbrack}}} = 0.} & (6.9) \end{matrix}$

Let us assume

(1) the functions ψ_(j)(y) and their derivatives are uniformly bounded for y≧0,

(2) ψ₀(y)=0 for y≧0. We can translate these assumptions into discrete transparent boundary conditions that the functions ψ_(j+1)(y)(j=0,1, . . . ) satisfy at y=0. A newly created algorithm computes these conditions recursively. First, define $\begin{matrix} {{p_{j + 1} = {{{e^{{- i}\quad{\pi/4}}\left( \frac{8}{\left( {x_{j + 1} - x_{j}} \right)\left\lbrack {{k^{- 1}\left( x_{j} \right)} + {k^{- 1}\left( x_{j + 1} \right)}} \right\rbrack} \right)}^{\frac{1}{2}}\quad{for}\quad j} = 0}},1,\ldots} & (6.10) \end{matrix}$ Second, define $\begin{matrix} {{g_{k}^{j} = {\frac{p_{k} - p_{j}}{p_{k} + p_{j}}\quad{for}\quad j}},{k = 1},2,\ldots} & (6.11) \end{matrix}$ Note that g_(k) ^(j) is real and |g_(k) ^(j) |<1 for all j, k.

If j=1, then we let $\begin{matrix} {a_{1}^{1} = {{\frac{1}{2}p_{1}^{- 2}\frac{\partial}{\partial y}{\psi_{1}(0)}\quad{and}\quad a_{2}^{1}} = {{- \frac{1}{2}}p_{1}^{- 2}\frac{\partial}{\partial y}{{\psi_{1}(0)}.}}}} & (6.12) \end{matrix}$ If j>1, then we compute the coefficients a₁ ^(j), . . . , a_(j+1) ^(j) from the coefficients a₁ ^(j−1), . . . , a_(j) ^(j−1) and from the derivatives σψ_(j−1)(ο)/σy and σψ_(j)(ο)/σy in two steps.

(1) Compute ā₁ ^(j−1), . . . , ā_(j) ^(j−1) and {overscore (b)}₁ ^(j−1), . . . , {overscore (b)}_(j) ^(j−1) from the following recurrence formulas. ā₁ ^(j−1)=a₁ ^(j−1),  (6.13) ā_(k) ^(j−1)=ā_(k) ^(j−1)−g_(k−1) ^(j)ā_(k−1) ^(j−1) for k=2, . . . ,j,  (6.14) {overscore (b)}₁ ^(j−1)=g₁ ^(j)ā₁ ^(j−1),  (6.15) {overscore (b)}_(k) ^(j−1)=g_(k) ^(j)ā_(k) ^(j−1)+ā_(k−1) ^(j−1)−g_(k−1) ^(j){overscore (b)}_(k−1) ^(j−1) for k=2, . . . ,j.  (6.16)

(2). Compute the coefficients a₁ ^(j), . . . , a_(j+1) ^(j) as follows. $\begin{matrix} {{a_{1}^{j} = {{{- \frac{1}{2}}a_{1}^{j - 1}} - {\frac{1}{2}g_{1}^{j}{\overset{\_}{b}}_{1}^{j - 1}}}},} & (6.17) \\ {{a_{k}^{j} = {{{{- \frac{1}{2}}a_{k}^{j - 1}} - {\frac{1}{2}g_{k}^{j}{\overset{\_}{b}}_{k}^{j - 1}} - {\frac{1}{2}{\overset{\_}{b}}_{k - 1}^{j - 1}\quad{for}\quad k}} = 2}},\ldots\quad,{j - 1},} & (6.18) \end{matrix}$

(Note: omit the preceding equation if j=2.) $\begin{matrix} {{a_{j}^{j} = {{\frac{1}{2}{p_{j}^{- 2}\left\lbrack {{\frac{\partial\quad}{\partial y}{\psi_{j - 1}(0)}} + {\frac{\partial\quad}{\partial y}{\psi_{j}(0)}}} \right\rbrack}} - {\frac{1}{2}a_{j}^{j - 1}} - {\frac{1}{2}a_{j}^{j - 1}} - {\frac{1}{2}{\overset{\_}{b}}_{j - 1}^{j - 1}}}},} & (6.19) \\ {a_{j + 1}^{j} = {{\frac{1}{2}{p_{j}^{- 2}\left\lbrack {{\frac{\partial\quad}{\partial y}{\psi_{j - 1}(0)}} + {\frac{\partial\quad}{\partial y}{\psi_{j}(0)}}} \right\rbrack}} - {\frac{1}{2}{{\overset{\_}{b}}_{j}^{j - 1}.}}}} & (6.20) \end{matrix}$ The coefficients ā_(j+1) ^(j) are used to compute the discrete transparent boundary conditions as follows. $\begin{matrix} {{{{\frac{\partial\quad}{\partial y}{\psi_{j}(0)}} + {p_{j + 1}{\psi_{j}(0)}}} = 0},} & (6.21) \\ {{{{\frac{\partial\quad}{\partial y}{\psi_{j + 1}(0)}} + {p_{j + 1}{\psi_{j + 1}(0)}}} = {{{{- \frac{\partial\quad}{\partial y}}{\psi_{j}(0)}} - {p_{j + 1}{\psi_{j}(0)}} + {2p_{j + 1}^{2}{\overset{\_}{a}}_{j + 1}^{j}\quad{for}\quad j}} = 1}},2,\ldots} & (6.22) \end{matrix}$

This algorithm can be greatly simplified if the quantities p_(j) are independent of the index j. This happens when the range increments x_(j+1)−x_(j) and the reference wavenumber values k(x_(j)) are independent of the index j. Let us assume that x_(j)=jΔx and k(x_(j))=k₀ for j=0,1, . . . , where the increment Δx and the wavenumber k₀ are fixed. Then it follows from eq. (6.10) that p_(j)=p₀ for every j, where p ₀ =e ^(ix/4)2(k ₀ /Δx)^(1/2).  (6.23) In addition, eq.(6.11) implies that g_(k) ^(j)=0 for all j,k. Thus, eqs.(6.13,6.14) become ā_(k) ^(j−1)=a_(k) ^(j−1) for k=1, . . . ,j,  (6.24) eq.(6.15) becomes {overscore (b)}₁ ^(j−1)=0 and eq.(6.16) becomes {overscore (b)}_(k) ^(j−1)=ā_(k−1) ^(j−1) for k=2, . . . ,j.  (6.25)

We can use these equations to eliminate the coefficients ā_(k) ^(j−1) and {overscore (b)}_(k) ^(j−1) from eqs.(6.17-6.20) if we extend the definition of the coefficients a_(k) ^(j) as follows. Let a⁻¹ ^(j)=0 for j=0,1, . . . ,  (6.26) a₀ ^(j)=0 for j=0,1 . . .  (6.27) Using these extensions, we get the following simplified forms of eqs.(6.17-6.20). For j>1, $\begin{matrix} {{a_{k}^{j} = {{{- \frac{1}{2}}a_{k}^{j - 1}} - {\frac{1}{2}a_{k}^{j - 1}} - {\frac{1}{2}a_{k - 2}^{j - 1}}}},{{{for}\quad k} = 1},\ldots\quad,{j - 1},} & (6.28) \\ {{a_{j}^{j} = {{\frac{1}{2}{p_{0}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{\psi_{j - 1}(0)}} + {\frac{\partial}{\partial y}{\psi_{j}(0)}}} \right\rbrack}} - a_{j}^{j - 1} - {\frac{1}{2}a_{j - 2}^{j - 1}}}},} & (6.29) \\ {a_{j + 1}^{j} = {{{- \frac{1}{2}}{p_{0}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{\psi_{j - 1}(0)}} + {\frac{\partial}{\partial y}{\psi_{j}(0)}}} \right\rbrack}} - {\frac{1}{2}{a_{j - 1}^{j - 1}.}}}} & (6.30) \end{matrix}$

We also define a₁ ⁰=0. This extension has the following useful consequences.

(1) We can let j=1 in eqs.(6.29,6.30), which become eq.(6.12) because σψ₀(0)/σy=0.

(2) Because of eq.(6.24), we can get eq.(6.21) from eq.(6.22) by letting j=0. Therefore, we can write these discrete transparent boundary conditions as follows. $\begin{matrix} {{{{\frac{\partial}{\partial y}{\psi_{j + 1}(0)}} + {p_{0}{\psi_{j + 1}(0)}}} = {{{- \frac{\partial}{\partial y}}{\psi_{j}(0)}} - {p_{0}{\psi_{j}(0)}} + {2p_{0}^{2}a_{j + 1}^{j}}}}{{{{for}\quad j} = 0},1,\ldots}} & (6.31) \end{matrix}$

(3) An alternative form of these conditions is obtained by iterating eq.(6.31) backwards. The initial conditions ψ₀(0)=0 and σψ₀(0)/σy=0 are used to stop the iteration. $\begin{matrix} {{{{\frac{\partial}{\partial y}{\psi_{j + 1}(0)}} + {p_{0}{\psi_{j + 1}(0)}}} = {2p_{0}^{2}{\sum\limits_{k = 0}^{j}{\left( {- 1} \right)^{j - k}a_{k + 1}^{k}}}}}\quad{{{{for}\quad j} = 0},1,\ldots}} & (6.32) \end{matrix}$

When eqs. (6.286.30) are evaluated in finite-precision computer arithmetic, roundoff errors in the additions tend to accumulate. However, we can use compensated summation to identify these errors and to cancel them to working precision. Let us describe how this is done in a simple case.

Consider the following sequence of differential boundary values σψ_(j)(0)/σy. $\begin{matrix} {{\frac{\partial}{\partial y}{\psi_{j}(0)}} = \left\{ \begin{matrix} {2p_{0}^{2}} & {{{{if}\quad j} = 1},} \\ 0 & {{{if}\quad j} \neq 1.} \end{matrix} \right.} & (6.33) \end{matrix}$ Since a₁ ⁰=0, it follows from eq. (6.32) that ψ₁(0)=−2p₀ and $\begin{matrix} {{{\psi_{j + 1}(0)} = {{2p_{0}{\sum\limits_{k = 1}^{j}{\left( {- 1} \right)^{j - k}a_{k + 1}^{k}\quad{for}\quad j}}} = 1}},2,\ldots} & (6.34) \end{matrix}$ Thus, ψ_(j+1)(0)=−2p₀σ_(j) where $\begin{matrix} {{\sigma_{j} = {{- {\sum\limits_{k = 1}^{j}{\left( {- 1} \right)^{j - k}a_{k + 1}^{k}\quad{for}\quad j}}} = 1}},2,\ldots\quad,} & (6.36) \end{matrix}$ which can be calculated recursively as follows. σ_(j)=−a_(j+1) ^(j)−σ_(j−1) for j=2,3, . . .  (6.36)

We can compute the coefficients a_(k) ^(j) for j=1, 2, 3 exactly from eqs.(6.28-6.30). Using the input values (6.33), we find that $\begin{matrix} \left. \begin{matrix} \begin{matrix} {a_{1}^{1} = 1} & {{a_{2}^{1} = {- 1}},} \end{matrix} \\ \begin{matrix} {a_{1}^{2} = {- \frac{1}{2}}} & {a_{2}^{2} = 2} & {{a_{3}^{2} = {- \frac{3}{2}}},} \end{matrix} \\ \begin{matrix} {a_{1}^{3} = \frac{1}{4}} & {a_{2}^{3} = {- 1}} & {a_{3}^{3} = \frac{7}{4}} & {a_{4}^{3} = {- 1.}} \end{matrix} \end{matrix} \right\} & (6.37) \end{matrix}$ Next, using the input values (6.33), we write eqs.(6.28-6.30) for j≧4 as follows. $\begin{matrix} {{a_{1}^{j} = {{- \frac{1}{2}}a_{1}^{j - 1}}},} & (6.38) \\ {{a_{2}^{j} = {{- \frac{1}{2}}a_{2}^{j - 1}}},} & (6.39) \\ {{a_{k}^{j} = {{{{- \frac{1}{2}}a_{k}^{j - 1}} - {\frac{1}{2}a_{k - 2}^{j - 1}\quad{for}\quad k}} = 3}},\ldots\quad,{j - 1},} & (6.40) \\ {{a_{j}^{j} = {{- a_{j}^{j - 1}} - {\frac{1}{2}a_{j - 2}^{j - 1}}}},} & (6.41) \\ {a_{j + 1}^{j} = {{- \frac{1}{2}}{a_{j - 1}^{j - 1}.}}} & (6.42) \end{matrix}$ We use the values of a₁ ³, a₂ ³, a₃ ³, a₄ ³ from eq.(6.37) to start these equations when j=4. We use eq.(6.36) to compute σ_(j). We start this equation with σ₃=½, which follows from eqs.(6.35, 6.37).

Compensated summation is based on the following principle. Let û and {circumflex over (υ)} denote real floating-point numbers, and let fl(û ±{circumflex over (υ)}) denote the floating-point sum and difference of û and {circumflex over (υ)}. Define the following floating-point numbers: ŝ=fl(û+{circumflex over (υ)}),  (6.43) û′=fl(ŝ−{circumflex over (υ)}),  (6.44) {circumflex over (υ)}″=fl(ŝ−{circumflex over (υ)}û′),  (6.45) {circumflex over (δ)}û′(û)=fl(û−û′),  (6.46) {circumflex over (δ)}″({circumflex over (υ)})=fl({circumflex over (υ)}−{circumflex over (υ)}″),  (6.47) Δ(û,{circumflex over (υ)})=fl({circumflex over (δ)}′(û)+{circumflex over (δ)}″({circumflex over (υ)})).  (6.48) The quantity Δ(û,{circumflex over (υ)}) is a very good estimate of the roundoff error(û+{circumflex over (υ)})−ŝ. Note that the arithmetic operations in eqs.(6.43-6.48) carry over directly to complex floating-point arithmetic because the operations on the real parts of the complex numbers and the operations on the imaginary parts of these numbers are identical but independent.

First, let us apply compensated summation to eq.(6.40). We want to compute $\begin{matrix} {{a_{k}^{j} = {{{- \left( {{\frac{1}{2}a_{k}^{j - 1}} + {\frac{1}{2}a_{k - 2}^{j - 1}}} \right)}\quad{for}\quad k} = 3}},\ldots\quad,{j - 1},} & (6.49) \end{matrix}$ but we actually compute $\begin{matrix} {{a_{k}^{j} = {{{- {{fl}\left( {{\frac{1}{2}a_{k}^{j - 1}} + {\frac{1}{2}a_{k - 2}^{j - 1}}} \right)}}\quad{for}\quad k} = 3}},\ldots\quad,{j - 1.}} & (6.50) \end{matrix}$ Let $\begin{matrix} {{\delta_{k}^{j} = {\left( {{\frac{1}{2}{\hat{a}}_{k}^{j - 1}} + {\frac{1}{2}{\hat{a}}_{k - 2}^{j - 1}}} \right) - {{fl}\left( {{\frac{1}{2}{\hat{a}}_{k}^{j - 1}} + {\frac{1}{2}{\hat{a}}_{k - 2}^{j - 1}}} \right)}}}\quad{{{{for}\quad k} = 3},\ldots\quad,{j - 1},}} & (6.51) \end{matrix}$ and ε_(k) ^(j)=a_(k) ^(j)−â_(k) ^(j) for k=1, . . . ,j+1.  (6.52) Now, eqs.(6.49-6.52) imply that $\begin{matrix} \begin{matrix} {a_{k}^{j} = {{- \left( {{\frac{1}{2}a_{k}^{j - 1}} + {\frac{1}{2}a_{k - 2}^{j - 1}}} \right)} + \delta_{k}^{j}}} \\ {= {{- \left( {{\frac{1}{2}a_{k}^{j - 1}} - {\frac{1}{2}\varepsilon_{k}^{j - 1}} + {\frac{1}{2}a_{k - 2}^{j - 1}} - {\frac{1}{2}\varepsilon_{k - 2}^{j - 1}}} \right)} + \delta_{k}^{j}}} \\ {= {{- \left( {{\frac{1}{2}a_{k}^{j - 1}} + {\frac{1}{2}a_{k - 2}^{j - 1}}} \right)} + {\frac{1}{2}\varepsilon_{k}^{j - 1}} + {\frac{1}{2}\varepsilon_{k - 2}^{j - 1}} + \delta_{k}^{j}}} \\ {{= {{a_{k}^{j} + {\frac{1}{2}\varepsilon_{k}^{j - 1}} + {\frac{1}{2}\varepsilon_{k - 2}^{j - 1}} + {\delta_{k}^{j}\quad{for}\quad k}} = 3}},\ldots\quad,{j - 1.}} \end{matrix} & (6.53) \end{matrix}$ Hence, $\begin{matrix} {{\varepsilon_{k}^{j} = {{\left( {{\frac{1}{2}\varepsilon_{k}^{j - 1}} + {\frac{1}{2}\varepsilon_{k - 2}^{j - 1}} + \delta_{k}^{j}} \right)\quad{for}\quad k} = 3}},\ldots\quad,{j - 1.}} & (6.54) \end{matrix}$ In practice we compute an approximation to ε_(k) ^(j) that we call {circumflex over (ε)}_(k) ^(j): $\begin{matrix} {{{\hat{\varepsilon}}_{k}^{j} = {{{- {{fl}\left( {{\frac{1}{2}{\hat{\varepsilon}}_{k}^{j - 1}} + {\frac{1}{2}{\hat{\varepsilon}}_{k - 2}^{j - 1}} + \Delta_{k}^{j}} \right)}}\quad{for}\quad k} = 3}},\ldots\quad,{j - 1},} & (6.55) \end{matrix}$ where {circumflex over (Δ)}_(k) ^(j) is the estimate for δ_(k) ^(j) that we get with $\hat{u} = {\frac{1}{2}{\hat{a}}_{k}^{j - 1}}$ and $\hat{v} = {\frac{1}{2}a_{k - 2}^{j - 1}}$ in eq. (6.48). Note that the RHS of eq.(6.55) is ambiguous because the associative law does not hold for floating-point addition. However, we shall ignore the small influence that the order of these additions may have on the numerical value of the total sum.

Next, we apply compensated summation to eq.(6.41). We want to compute $\begin{matrix} {{a_{j}^{j} = {- \left( {a_{j}^{j - 1} + {\frac{1}{2}a_{j - 2}^{j - 1}}} \right)}},} & (6.56) \end{matrix}$ but we actually compute $\begin{matrix} {{\hat{a}}_{j}^{j} = {{- f}\quad{{l\left( {{\hat{a}}_{j}^{j - 1} + {\frac{1}{2}{\hat{a}}_{j - 2}^{j - 1}}} \right)}.}}} & (6.57) \end{matrix}$ Let $\begin{matrix} {\delta_{j}^{j} = {\left( {{\hat{a}}_{j}^{j - 1} + {\frac{1}{2}{\hat{a}}_{j - 2}^{j - 1}}} \right) - {{{fl}\left( {{\hat{a}}_{j}^{j - 1} + {\frac{1}{2}{\hat{a}}_{j - 2}^{j - 1}}} \right)}.}}} & (6.58) \end{matrix}$ Now, eqs.(6.52,6.56-6.58) imply that $\begin{matrix} \begin{matrix} {{\hat{a}}_{j}^{j} = {\left( {{\hat{a}}_{j}^{j - 1} + {\frac{1}{2}{\hat{a}}_{j - 2}^{j - 1}}} \right) + \delta_{j}^{j}}} \\ {= {{- \left( {a_{j}^{j - 1} - \varepsilon_{j}^{j - 1} + {\frac{1}{2}a_{j - 2}^{j - 1}} - {\frac{1}{2}\varepsilon_{j - 2}^{j - 1}}} \right)} + \delta_{j}^{j}}} \\ {= {{- \left( {a_{j}^{j - 1} + {\frac{1}{2}a_{j - 2}^{j - 1}}} \right)} + \varepsilon_{j}^{j - 1} + {\frac{1}{2}\varepsilon_{j - 2}^{j - 1}} + \delta_{j}^{j}}} \\ {= {a_{j}^{j} + \varepsilon_{j}^{j - 1} + {\frac{1}{2}\varepsilon_{j - 2}^{j - 1}} + {\delta_{j}^{j}.}}} \end{matrix} & (6.59) \end{matrix}$ Hence, $\begin{matrix} {\varepsilon_{j}^{j} = {- {\left( {\varepsilon_{j}^{j - 1} + {\frac{1}{2}\varepsilon_{j - 2}^{j - 1}} + \delta_{j}^{j}} \right).}}} & (6.60) \end{matrix}$ In practice we compute an approximation to ε_(j) ^(j) that we call {circumflex over (ε)}_(j) ^(j): $\begin{matrix} {{{\hat{\varepsilon}}_{j}^{j} = {- {{fl}\left( {{\hat{\varepsilon}}_{j}^{j - 1} + {\frac{1}{2}{\hat{\varepsilon}}_{j - 2}^{j - 1}} + {\hat{\Delta}}_{j}^{j}} \right)}}},} & (6.61) \end{matrix}$ where {circumflex over (Δ)}_(j) ^(j) is the estimate for δ_(j) ^(j) that we get with û=â_(j) ^(j−1) and $\hat{v} = {\frac{1}{2}a_{j - 2}^{j - 1}}$ in eq.(6.48). The remarks that follow eq.(6.55) apply to eq. (6.61) as well.

The treatment of eqs.(6.38,6.39,6.42) is simple because they involve no summations. Since we want to compute ${a_{1}^{j} = {{- \frac{1}{2}}a_{1}^{j - 1}}},{a_{2}^{j} = {{- \frac{1}{2}}a_{2}^{j - 1}}},{a_{j + 1}^{j} = {{- \frac{1}{2}}a_{j - 1}^{j - 1}}},$ the computed values of these coefficients satisfy the same equations: $\begin{matrix} {{{\hat{a}}_{1}^{j} = {{- \frac{1}{2}}{\hat{a}}_{1}^{j - 1}}},{{\hat{a}}_{2}^{j} = {{- \frac{1}{2}}{\hat{a}}_{2}^{j - 1}}},{{\hat{a}}_{j + 1}^{j} = {{- \frac{1}{2}}{{\hat{a}}_{j - 1}^{j - 1}.}}}} & (6.62) \end{matrix}$ Now the errors ε₁ ^(j), ε₂ ^(j), ε_(j+1) ^(j) also satisfy these equations: $\begin{matrix} {{\varepsilon_{1}^{j} = {{- \frac{1}{2}}\varepsilon_{1}^{j - 1}}},{\varepsilon_{2}^{j} = {{- \frac{1}{2}}\varepsilon_{2}^{j - 1}}},{\varepsilon_{j + 1}^{j} = {{- \frac{1}{2}}{\varepsilon_{j - 1}^{j - 1}.}}}} & (6.63) \end{matrix}$ Hence, we use the same equations to compute the approximations {circumflex over (ε)}₁ ^(j), {circumflex over (ε)}₂ ^(j), {circumflex over (ε)}_(j+1) ^(j) to these errors: $\begin{matrix} {{{\hat{\varepsilon}}_{1}^{j} = {{- \frac{1}{2}}{\hat{\varepsilon}}_{1}^{j - 1}}},{{\hat{\varepsilon}}_{2}^{j} = {{- \frac{1}{2}}{\hat{\varepsilon}}_{2}^{j - 1}}},{{\hat{\varepsilon}}_{j + 1}^{j} = {{- \frac{1}{2}}{{\hat{\varepsilon}}_{j - 1}^{j - 1}.}}}} & (6.64) \end{matrix}$

The approximate corrections {circumflex over (ε)}_(k) ³ vanish because the computed coefficients â_(k) ³ are exact. We use the values {circumflex over (ε)}₁ ³=0, {circumflex over (ε)}₂ ³=0, {circumflex over (ε)}₃ ³=0, {circumflex over (ε)}₄ ³=0 to start eqs.(6.55,6.61,6.64) when j=4. Since {circumflex over (ε)}₁ ³=0 and {circumflex over (ε)}₂ ³=0, the first two parts of eq.(6.64) imply that {circumflex over (ε)}₁ ^(j)=0 and {circumflex over (ε)}₂ ^(j)=0 for all j≧4. Therefore, we need only the last part of eq.(6.64).

Finally, we consider {circumflex over (σ)}_(j) and develop an approximate correction to it. We want to compute σ_(j)=−(a_(j+1) ^(j)+σ_(j−1)) for j≧4,  (6.65) with σ₃=½. In practice, we compute {circumflex over (σ)}_(j)=−fl(â_(j+1) ^(j)+{circumflex over (σ)}_(j−1)) for j≧4,  (6.66) with {circumflex over (σ)}₃=½. Define e_(j)=σ_(j)−{circumflex over (σ)}_(j) for j≧3.  (6.67) By the same approach that lead to eqs.(6.54,6.60) we find that e_(j)=−(ε_(j+1) ^(j)+e_(j−1)+δ_(j)) for j≧4,  (6.68) where δ_(j)=â_(j+1) ^(j)+{circumflex over (σ)}_(j−1)−fl(â_(j+1) ^(j)+{circumflex over (σ)}_(j−1)).  (6.69) Although we could use eq.(6.48) to estimate δ_(j) we shall neglect it against ε_(j+1) ^(j). Therefore, we compute an approximate correction ê_(j) to {circumflex over (σ)}_(j) by the following recurrence formula. {circumflex over (ε)}_(j)=−fl({circumflex over (68 )}_(j+1) ^(j)+ê_(j−1)) for j≧4,  (6.70) where ê₃=0.

Let us summarize the computational forms of eqs.(6.36,6.38-6.42) for j≧4 that use compensated summation to identify roundoff errors and to cancel them to working precision. $\begin{matrix} \left. \begin{matrix} {{{\hat{a}}_{1}^{j} = {{- \frac{1}{2}}{\hat{a}}_{1}^{j - 1}}},} \\ {{{\hat{a}}_{2}^{j} = {{- \frac{1}{2}}{\hat{a}}_{1}^{j - 1}}},} \\ {{{{for}\quad k} = 3},\ldots\quad,{{j - 1}:}} \\ {{{\hat{a}}_{k}^{j} = {{{- {{fl}\left( {{\frac{1}{2}{\hat{a}}_{k}^{j - 1}} + {\frac{1}{2}{\hat{a}}_{k - 2}^{j - 1}}} \right)}}\quad{and}\quad{\hat{\varepsilon}}_{k}^{- j}} = {- {{fl}\left( {{\frac{1}{2}{\hat{\varepsilon}}_{k}^{j - 1}} + {\frac{1}{2}{\hat{\varepsilon}}_{k - 2}^{j - 1}} + {\hat{\Delta}}_{k}^{j}} \right)}}}},} \\ {{{\hat{a}}_{j}^{j} = {{{- {{fl}\left( {{\hat{a}}_{j}^{j - 1} + {\frac{1}{2}{\hat{a}}_{j - 2}^{j - 1}}} \right)}}\quad{and}\quad{\hat{\varepsilon}}_{j}^{- j}} = {- {{fl}\left( {{\hat{\varepsilon}}_{j}^{j - 1} + {\frac{1}{2}{\hat{\varepsilon}}_{j - 2}^{j - 1}} + {\hat{\Delta}}_{j}^{j}} \right)}}}},} \\ {{{\hat{a}}_{j + 1}^{j} = {{\frac{1}{2}{\hat{a}}_{j - 1}^{j - 1}\quad{and}\quad{\hat{\varepsilon}}_{j + 1}^{- j}} = {{- \frac{1}{2}}{\hat{\varepsilon}}_{j - 1}^{j - 1}}}},} \\ {{\hat{\sigma}}_{j} = {{{- {{fl}\left( {{\hat{a}}_{j + 1}^{j} + {\hat{\sigma}}_{j - 1}} \right)}}\quad{and}\quad{\hat{e}}_{j}} = {- {{{fl}\left( {{\hat{\varepsilon}}_{j + 1}^{- j} + {\hat{e}}_{j - 1}} \right)}.}}}} \end{matrix} \right\} & (6.71) \end{matrix}$

We approximate σ_(j) by the sum {circumflex over (σ)}_(j)+ê_(j).

The number of arithmetic operations in eq. (6.28) of the simple algorithm is proportional to j−1, which increases with range. We can limit the size of this operation count by truncating the algorithm. First, we choose an integer L≧3, which will be called the length of the truncated algorithm. If j≧L, then we discard all the coefficients a_(k) ^(j−1) such that k≦j−L. We effectively equate these coefficients to 0, which lets us start the recursion in eq.(6.28) for the truncated algorithm at k=j−L+1. Since the coefficient a_(j−L+1) ^(j) is not used in subsequent steps, however, there is no need to evaluate it. For this reason, if j≧L, then we can start the recursion in eq.(6.28) for the truncated algorithm at k=j−L+2. Now the operation count is proportional to L−2, which is fixed.

We shall state the truncated algorithm in terms of new coefficients a_(l) ^(j−1) that are defined for j=1,2, . . . and l=1, . . . ,L as follows.

(1) Let a_(l) ⁰=0 for l=1, . . . ,L.  (6.72)

(2) Let a₀ ^(j)=0 for j=0,1, . . .  (6.73)

(3) For j=1,2, . . . , $\begin{matrix} {{\alpha_{l}^{j} = {{{- \frac{1}{2}}\alpha_{l + 1}^{j - 1}} - {\frac{1}{2}\alpha_{l - 1}^{j - 1}}}},{{{for}\quad l} = 1},\ldots\quad,{L - 2},} & (6.74) \\ {{\alpha_{L - 1}^{j} = {{\frac{1}{2}{p_{0}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{\psi_{j - 1}(0)}} + {\frac{\partial}{\partial y}{\psi_{j}(0)}}} \right\rbrack}} - \alpha_{L}^{j - 1} - {\frac{1}{2}\alpha_{L - 2}^{j - 1}}}},} & (6.75) \\ {\alpha_{L}^{j} = {{\frac{1}{2}{p_{0}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{\psi_{j - 1}(0)}} + {\frac{\partial}{\partial y}{\psi_{j}(0)}}} \right\rbrack}} - {\frac{1}{2}{\alpha_{L - 1}^{j - 1}.}}}} & (6.76) \end{matrix}$ The transparent boundary conditions associated with this algorithm are $\begin{matrix} {{{{\frac{\partial}{\partial y}{\psi_{j + 1}(0)}} + {p_{0}{\psi_{j + 1}^{(L)}(0)}}} = {{{{- \frac{\partial}{\partial y}}{\psi_{j}(0)}} - {p_{0}{\psi_{j}^{(L)}(0)}} + {2p_{0}^{2}\alpha_{L}^{j}\quad{for}\quad j}} = 0}},1,\ldots} & (6.77) \end{matrix}$ The superscript “L” indicates that the sequence of boundary values {ψ_(j) ^((L))(0)}_(j=0) ^(∞) is produced from the sequence of boundary derivatives {σψ_(j)(0)/σy}_(j=0) ^(∞) by the truncated algorithm of length L. These boundary values are only approximate.

An alternative form of these conditions is obtained by iterating eq.(6.77) backwards. The initial conditions ψ₀ ^((L))(0)=0 and σψ₀(0)/σy=0 are used to stop the iteration. The result is like eq.(6.32). $\begin{matrix} {{{{\frac{\partial}{\partial y}{\psi_{j + 1}(0)}} + {p_{0}{\psi_{j + 1}^{(L)}(0)}}} = {{2p_{0}^{2}{\sum\limits_{k = 0}^{j}{\left( {- 1} \right)^{j - k}\alpha_{L}^{k}\quad{for}\quad j}}} = 0}},1,\ldots} & (6.78) \end{matrix}$

Suppose the sequence of differential boundary values σψ_(j)(0)/σy in eq.(6.33) is used as input to the truncated algorithm. The errors in the approximate boundary values ψ_(j) ^((L))(0) relative to the exact boundary values ψ_(j)(0) can be characterized as follows. Fix ξ>0, and let {j₁, j₂, . . . j_(l), . . . } be a sequence of integers such that j_(l)˜ξl² as l→∞. Then the relative error in ψ_(jL) ^((L))(0) converges to a finite limit as L→∞ that is exponentially small if ξ is small. $\begin{matrix} {0 < {\lim\limits_{L\rightarrow\infty}\left( {1 - {{\psi_{jL}^{(L)}(0)}/{\psi_{jL}(0)}}} \right)} < {2{c^{{- 2}/ɛ}.}}} & (6.79) \end{matrix}$

Suppose the index j and the length L are given, and let ξ=jL⁻². A practical way to interpret eq.(6.79) is that 1−ψ_(j) ^((L))(0)/ψ_(j)(0)≈2 exp(−2/ξ) if L >>1.

Let u denote the unit roundoff in a finite-precision floating-point computer arithmetic. Two real numbers a and b are represented by the same floating-point number if |1−a/b|<u. These numbers are equivalent in this computer arithmetic. If L>[ln(2/u)j/2]^(1/2), then 2 exp(−2/ξ)<u. Hence, ψ_(j) ^((L))(0) and ψ_(j)(0) are equivalent in the computer arithmetic if L>>1.

Now let J be an integer, and let L_(J) be the smallest integer greater than [ln(2/u)J/2]^(1/2). L_(j)=[|ln(2/u)J/2)^(1/2)]=ο(J^(1/2)).  (6.80) If L_(J)>>1, then ψ_(j) ^((L) ^(J) ⁾(0) and ψ_(j)(0) are equivalent in the computer arithmetic for all j≦J. The number of arithmetic operations to compute each ψ_(j) ^((L) ^(J) ⁾(0) is proportional to L_(J), but the number of arithmetic operations to compute each ψ^(j)(0) is proportional to j. The ratio ρ(J) of the total operation count for computing the subsequence {ψ_(j) ^((L) ^(J) ⁾(0)}_(j=1) ^(J) to the total operation count for computing the subsequence {ψ_(j)(0)}_(j=1) ^(J) is approximately $\begin{matrix} {{{\rho(J)} \approx {J\quad{L_{J}/\frac{1}{2}}J^{2}}} = {{2{L_{J}/J}} = {{O\left( J^{{- 1}/2} \right)}.}}} & (6.81) \end{matrix}$

Therefore, if L_(J)>>1, then the truncated algorithm of length L_(J) is more efficient for computing the first J boundary values to working precision than the complete algorithm. As an example of practical importance, consider IEEE double-precision floating-point arithmetic. Since u=2⁻⁵³ in this arithmetic, L_(J)≈4.326J^(1/2) for IEEE double-precision floating-point arithmetic.  (6.82) These criteria are valid for all sequences of differential boundary values σψ_(j)(0)/σy. However, let us continue to use the sequence in eq.(6.33) to explain how to program the truncated algorithm with compensated summation.

Assume that L>4. We can calculate the coefficients a_(l) ^(j) for j=1,2,3 exactly. $\begin{matrix} \left. \begin{matrix} \begin{matrix} {\alpha_{l}^{1} = {{0\quad{for}\quad l} \leq {L - 2}}} & {\alpha_{L - 1}^{1} = 1} & {{\alpha_{L}^{1} = {- 1}},} \end{matrix} \\ \begin{matrix} {\alpha_{l}^{2} = {{0\quad{for}\quad l} \leq {L - 3}}} & {\alpha_{L - 2}^{2} = {- \frac{1}{2}}} & {\alpha_{L - 1}^{2} = 2} & {{\alpha_{L}^{2} = {- \frac{3}{2}}},} \end{matrix} \\ \begin{matrix} {\alpha_{l}^{3} = {{0\quad{for}\quad l} \leq {L - 4}}} & {\alpha_{L - 3}^{3} = \frac{1}{4}} & {\alpha_{L - 2}^{3} = {- 1}} & {\alpha_{L - 1}^{3} = {{\frac{7}{4}\quad\alpha_{L}^{3}} = {- 1.}}} \end{matrix} \end{matrix} \right\} & (6.83) \end{matrix}$ Next, we use eqs.(6.33,6.73) to write eqs.(6.74-6.76) for j≧4 as follows. $\begin{matrix} {{\alpha_{1}^{j} = {{- \frac{1}{2}}\alpha_{2}^{j - 1}}},} & (6.84) \\ {{\alpha_{l}^{j} = {{{{- \frac{1}{2}}\alpha_{l + 1}^{j - 1}} - {\frac{1}{2}\alpha_{l - 1}^{j - 1}\quad{for}\quad l}} = 2}},\ldots\quad,{L - 2},} & (6.85) \\ {{\alpha_{L - 1}^{j} = {{- \alpha_{L}^{j - 1}} - {\frac{1}{2}\alpha_{L - 2}^{j - 1}}}},} & (6.86) \\ {\alpha_{L}^{j} = {{- \frac{1}{2}}{\alpha_{L - 1}^{j - 1}.}}} & (6.87) \end{matrix}$

Since a_(L) ⁰=0, it follows from eq.(6.78) that ψ₁ ^((L))(0)=−2p₀ and $\begin{matrix} {{{\psi_{j + 1}^{(L)}(0)} = {{2\quad p_{0}{\sum\limits_{k = 1}^{j}{\left( {- 1} \right)^{j - k}\alpha_{L}^{k}\quad{for}\quad j}}} = 1}},2,{\ldots\quad.}} & (6.88) \end{matrix}$ Thus, ψ_(j+1) ^((L))(0)=−2p₀σ_(j) ^((L)) where $\begin{matrix} {{\sigma_{j}^{(L)} = {{- {\sum\limits_{k = 1}^{j}{\left( {- 1} \right)^{j - k}\alpha_{L}^{k}\quad{for}\quad j}}} = 1}},2,\ldots\quad,} & (6.89) \end{matrix}$ which can be calculated recursively as follows. σ_(j) ^((L))=a_(L) ^(j)−σ_(j−1) ^((L)) for j=2,3, . . .  (6.90) We use the values of a_(l) ³ for l=1, . . . ,L from eq.(6.83) to start eqs.(6.84-6.87) when j=4. We compute σ_(j) ^((L)) from eq.(6.90). We start this equation with σ₃ ^((L))={fraction (1/2)}, which follows from eqs.(6.83,6.89).

If eqs.(6.84-6.87,6.90) are programmed directly, then roundoff errors tend to accumulate. We can use compensated summation to counteract these roundoff errors in the same way that we use it to counteract roundoff errors for the complete algorithm. Let us denote the computed values of a_(k) ^(j) and σ_(j) ^((L)) by â_(k) ^(j) and {circumflex over (σ)}_(j) ^((L)), respectively. We list the corrected computational forms of eqs.(6.84-6.87,6.90) next. If j≧4, then $\left. \begin{matrix} {{{\hat{\alpha}}_{1}^{j} = {{{- \frac{1}{2}}{\hat{\alpha}}_{2}^{j - 1}\quad{and}\quad{\hat{\varepsilon}}_{1}^{j}} = {{- \frac{1}{2}}{\hat{\varepsilon}}_{2}^{j - 1}}}},} \\ {{{{for}\quad l} = 2},\ldots\quad,{L - {2\text{:}}}} \\ {{{\hat{\alpha}}_{l}^{j} = {{{- {{fl}\left( {{\frac{1}{2}{\hat{\alpha}}_{l + 1}^{j - 1}} + {\frac{1}{2}{\hat{\alpha}}_{l - 1}^{j - 1}}} \right)}}\quad{and}\quad{\hat{\varepsilon}}_{l}^{j}} = {- {{fl}\left( {{\frac{1}{2}{\hat{\varepsilon}}_{l + 1}^{j - 1}} + {\frac{1}{2}{\hat{\varepsilon}}_{l - 1}^{j - 1}} + {\hat{\Delta}}_{l}^{j}} \right)}}}},} \\ {{{\hat{\alpha}}_{L - 1}^{j} = {{{- {{fl}\left( {{\hat{\alpha}}_{L}^{j - 1} + {\frac{1}{2}{\hat{\alpha}}_{L - 2}^{j - 1}}} \right)}}\quad{and}\quad{\hat{\varepsilon}}_{L - 1}^{j}} = {- {{fl}\left( {{\hat{\varepsilon}}_{L}^{j - 1} + {\frac{1}{2}{\hat{\varepsilon}}_{L - 2}^{j - 1}} + {\hat{\Delta}}_{L - 1}^{j}} \right)}}}},} \\ {{{\hat{\alpha}}_{L}^{j} = {{{- \frac{1}{2}}{\hat{\alpha}}_{L - 1}^{j - 1}\quad{and}\quad{\hat{\varepsilon}}_{L}^{j}} = {{- \frac{1}{2}}{\hat{\varepsilon}}_{L - 1}^{j - 1}}}},} \\ {{\hat{\sigma}}_{j}^{(L)} = {{{- {{fl}\left( {{\hat{\alpha}}_{L}^{j} + {\hat{\sigma}}_{j - 1}^{(L)}} \right)}}\quad{and}\quad{\hat{\varepsilon}}_{j}^{(L)}} = {- {{{fl}\left( {{\hat{\varepsilon}}_{L}^{j} + {\hat{\varepsilon}}_{j - 1}^{(L)}} \right)}.}}}} \end{matrix} \right\}\quad(6.91)$ Using eqs.(6.43-6.48),we compute {circumflex over (Δ)}_(l) ^(j) for l=2, . . . ,L−2 as {circumflex over (Δ)}(û,{circumflex over (υ)}) with $\hat{u} = {\frac{1}{2}{\hat{\alpha}}_{l + 1}^{j - 1}}$ and ${\hat{\upsilon} = {\frac{1}{2}{\hat{\alpha}}_{l - 1}^{j - 1}}},$ and we compute {circumflex over (Δ)}_(L−1) ^(j) as {circumflex over (Δ)}(û,{circumflex over (υ)}) with û=â_(L) ^(j−1) and $\hat{\upsilon} = {\frac{1}{2}{{\hat{\alpha}}_{L - 2}^{j - 1}.}}$

These equations are developed like the ones that are collected as eq.(6.71). We start them when j=4 with the following initial values: $\begin{matrix} \begin{matrix} {{\hat{\alpha}}_{l}^{3} = {{0\quad{for}\quad l} \leq {L - 4}}} \\ {{\hat{\alpha}}_{L - 3}^{3} = \frac{1}{4}} \\ {{\hat{\alpha}}_{L - 2}^{3} = {- 1}} \\ {{\hat{\alpha}}_{L - 1}^{3} = \frac{7}{4}} \\ {{\hat{\alpha}}_{L}^{3} = {- 1}} \\ {{{\hat{\sigma}}_{3}^{(L)} = \frac{1}{2}},} \\ {{\hat{\varepsilon}}_{l}^{3} = {{0\quad{for}\quad l} \leq L}} \\ {{\hat{e}}_{3}^{(L)} = 0.} \end{matrix} & (6.92) \end{matrix}$ We approximate σ_(j) ^((L)) by the sum {circumflex over (σ)}_(j) ^((L)), which we denote by {circumflex over (σ)}_(j) ^((L).) 7. A Chebyshev tau method for semi-discrete parabolic equations.

If x_(j)=jΔx and k(x_(j))=k₀ for j=0,1, . . . , where the increment Δx and the wavenumber k₀ are fixed, then eq.(6.9) takes the following simple form. $\begin{matrix} {{{\psi_{j + 1}(y)} - {\psi_{j}(y)}} = {{{\left( {i\quad\Delta\quad{x/4}k_{0}} \right)\left\lbrack {{\frac{\partial^{2}}{\partial y^{2}}{\psi_{j}(y)}} + {\frac{\partial^{2}}{\partial y^{2}}{\psi_{j + 1}(y)}}} \right\rbrack}\quad{for}\quad j} = {\quad{0,1,\ldots}}}} & (7.1) \end{matrix}$ We can use Chebyshev interpolation polynomials to integrate this equation numerically over an interval −W≦y≦W. Let T_(k)(ξ) denote the k-th Chebyshev polynomial (of the first kind). For k=0,1, . . . these polynomials are defined by the identity T_(k)(ξ)=cos(kθ), where ξ=cos θ.  (7.2) Let N be a positive integer. A Chebyshev interpolation polynomial of degree N over the interval −W≦y≦W is a sum of the form $\begin{matrix} {{f(y)} = {\sum\limits_{k = 0}^{N}\quad{f_{k}{{T_{k}\left( {y/W} \right)}.}}}} & (7.3) \end{matrix}$

We shall approximate each function ψ_(j)(y) by a polynomial ƒ_(j)(y) of degree N, which we represent uniquely as a Chebyshev interpolation polynomial. $\begin{matrix} {{f_{j}(y)} = {{{\sum\limits_{k = 0}^{N}\quad{f_{j,k}{T_{k}\left( {y/W} \right)}\quad{for}}} - W} \leq y \leq {W.}}} & (7.4) \end{matrix}$ The first and second derivatives of ƒ_(j)(y) with respect to y are polynomials of degrees N−1 and N−2. We also represent these polynomials uniquely as Chebyshev interpolation polynomials. $\begin{matrix} {{{\frac{\partial}{\partial y}{f_{j}(y)}} = {{{\frac{1}{W}{\sum\limits_{k = 0}^{N - 1}\quad{f_{j,k}^{(1)}{T_{k}\left( {y/W} \right)}\quad{for}}}} - W} \leq y \leq W}},} & (7.5) \end{matrix}$ $\begin{matrix} {{\frac{\partial^{2}}{\partial y^{2}}{f_{j}(y)}} = {{{\frac{1}{W^{2}}{\sum\limits_{k = 0}^{N - 2}\quad{f_{j,k}^{(2)}{T_{k}\left( {y/W} \right)}\quad{for}}}} - W} \leq y \leq {W.}}} & (7.6) \end{matrix}$

Assume that the coefficients in these equations satisfy the following Chebyshev tau equations. $\begin{matrix} {{f_{{j + 1},k} - f_{j,k}} = {\left( {i\quad\Delta\quad{x/4}k_{0}W^{2}} \right)\left( {f_{{j + 1},k}^{(2)} + f_{j,k}^{(2)}} \right)\quad{for}\left\{ \begin{matrix} {{j = 0},1,{\ldots\quad{and}}} \\ {{k = 0},1,\ldots\quad,{N - 2.}} \end{matrix} \right.}} & (7.7) \end{matrix}$ Under this assumption, the polynomials ƒ_(j)(y) satisfy equations that are similar to eq.(7.1). For j=0,1, . . . and −W≦y≦W, $\begin{matrix} {{{{f_{j + 1}(y)} - {f_{j}(y)}} = {{\left( {i\quad\Delta\quad{x/4}k_{0}} \right)\left\lbrack {{\frac{\partial^{2}}{\partial y^{2}}{f_{j + 1}(y)}} + {\frac{\partial^{2}}{\partial y^{2}}{f_{j}(y)}}} \right\rbrack} + {e_{{j + 1},{N - 1}}{T_{N - 1}\left( {y/W} \right)}} + {e_{{j + 1},N}{T_{N}\left( {y/W} \right)}}}},} & (7.8) \end{matrix}$ where e_(j+1,N−1)=ƒ_(j+1,N−1) and e_(j+1,N)=ƒ_(j+1,N)−ƒ_(j,N).  (7.9)

We also define the following dimensionless and real-valued parameter. q=4k₀W²/Δz.  (7.10) Now the Chebyshev tau equations can be written in the following equivalent form. $\begin{matrix} {{f_{{j + 1},k}^{(2)} + {{iq}\quad f_{{j + 1},k}}} = {g_{j,k}\quad{for}\left\{ \begin{matrix} {{j = 0},1,{\ldots\quad{and}}} \\ {{k = 0},1,\ldots\quad,{N - 2},} \end{matrix} \right.}} & (7.11) \end{matrix}$ where g_(j,k)=−ƒ_(j,k) ⁽²⁾+iqƒ_(j,k).  (7.12) The quantities g_(j,k) can be computed recursively without the coefficients ƒ_(j+1,k) ⁽²⁾. $\begin{matrix} \begin{matrix} {g_{{j + 1},k} = {{- f_{{j + 1},k}^{(2)}} + {{iq}\quad f_{{j + 1},k}}}} \\ {= {{- \left( {f_{{j + 1},k}^{(2)} + {{iq}\quad f_{{j + 1},k}}} \right)} + {2{iq}\quad f_{{j + 1},k}}}} \\ {= {{- g_{j,k}} + {2{iq}\quad{f_{{j + 1},k}.}}}} \end{matrix} & (7.13) \end{matrix}$

It is possible also to eliminate the coefficients ƒ_(j+1,k) ⁽²⁾ from the Chebyshev tau equations. This leads to an alternative formulation of these equations in terms of the coefficients ƒ_(j+1,k) alone. The coefficients ƒ_(j+1,k) and ƒ_(j+1,k) ⁽²⁾ in eqs.(7.4,7.6) satisfy eq.(7.7) if and only if. $\begin{matrix} {{{f_{{j + 1},N} + {\frac{1}{4}\frac{1}{N\left( {N - 1} \right)}{iq}\quad f_{{j + 1},{N - 2}}}} = h_{j,N}},} & (7.14) \\ {{{f_{{j + 1},{N - 1}} + {\frac{1}{4}\frac{1}{\left( {N - 1} \right)\left( {N - 2} \right)}{iq}\quad f_{{j + 1},{N - 3}}}} = h_{j,{N - 1}}},} & (7.15) \\ {{{{\left( {1 - {\frac{1}{2}\frac{1}{\left( {N - 1} \right)\left( {N - 3} \right)}{iq}}} \right)f_{{j + 1},{N - 2}}} + {\frac{1}{4}\frac{1}{\left( {N - 2} \right)\left( {N - 3} \right)}{iq}\quad f_{{j + 1},{N - 4}}}} = h_{j,{N - 2}}},} & (7.16) \\ {{{{\left( {1 - {\frac{1}{2}\frac{1}{\left( {N - 2} \right)\left( {N - 4} \right)}{iq}}} \right)f_{{j + 1},{N - 3}}} + {\frac{1}{4}\frac{1}{\left( {N - 3} \right)\left( {N - 4} \right)}{iq}\quad f_{{j + 1},{N - 5}}}} = h_{j,{N - 3}}},} & (7.17) \\ {{{{\frac{1}{4}\frac{1}{\left( {k + 1} \right)k}{iq}\quad f_{{j + 1},{k + 2}}} + {\left( {1 - {\frac{1}{2}\frac{1}{k^{2} - 1}{iq}}} \right)f_{{j + 1},k}} + {\frac{1}{4}\frac{1}{k\left( {k - 1} \right)}{iq}\quad f_{{j + 1},{k - 2}}}} = {{h_{j,k}\quad{for}\quad k} = {N - 4}}},\ldots\quad,3,} & (7.18) \end{matrix}$  for k=N−4, . . . ,3,  (7.18) $\begin{matrix} {{{{\frac{1}{24}{iq}\quad f_{{j + 1},4}} + {\left( {1 - {\frac{1}{6}{iq}}} \right)f_{{j + 1},2}} + {\frac{1}{4}{iq}\quad f_{{j + 1},0}}} = h_{j,2}},} & (7.19) \end{matrix}$ for j=0,1, . . . , where $\begin{matrix} {{h_{j,N} = {\frac{1}{4}\frac{1}{N\left( {N - 1} \right)}g_{j,{N - 2}}}},} & (7.20) \\ {h_{j,{N - 1}} = {\frac{1}{4}\frac{1}{\left( {N - 1} \right)\left( {N - 2} \right)}g_{j,{N - 3},}}} & (7.21) \\ {h_{j,{N - 2}} = {{{- \frac{1}{2}}\frac{1}{\left( {N - 1} \right)\left( {N - 3} \right)}g_{j,{N - 2}}} + {\frac{1}{4}\frac{1}{\left( {N - 2} \right)\left( {N - 3} \right)}g_{j,{N - 4},}}}} & (7.22) \\ {h_{j,{N - 3}} = {{{- \frac{1}{2}}\frac{1}{\left( {N - 2} \right)\left( {N - 4} \right)}g_{j,{N - 3}}} + {\frac{1}{4}\frac{1}{\left( {N - 3} \right)\left( {N - 4} \right)}g_{j,{N - 5},}}}} & (7.23) \\ {{h_{j,k} = {{{\frac{1}{4}\frac{1}{\left( {k + 1} \right)k}g_{j,{k + 2}}} - {\frac{1}{2}\frac{1}{k^{2} - 1}g_{j,k}} + {\frac{1}{4}\frac{1}{k\left( {k - 1} \right)}g_{j,{k - 2}}\quad{for}\quad k}} = {N - 4}}},\ldots\quad,3,} & (7.24) \\ {h_{j,2} = {{\frac{1}{24}g_{j,4}} - {\frac{1}{6}g_{j,2}} + {\frac{1}{4}{g_{j,0}.}}}} & (7.25) \end{matrix}$ From this point on we assume that N is even. Note that

-   -   (1) there are N−1 equations in eqs.(7.14-7.19); ,     -   (2) the even-indexed coefficients ƒ_(j+1,N), ƒ_(j+1,N−2), . . .         ƒ_(j+1), ƒ_(j+1,0) are uncoupled from the odd-indexed         coefficients ƒ_(j+1,N−1), ƒ_(j+1,N−3), . . . ,ƒ_(j+1,1) in         eqs.(7.14-7.19);     -   (3) the coefficient f_(j+1,N) appears only in eq.(7.14), and the         coefficient ƒ_(j+1,N−1) appears only in eq.(7.15).

We adopt conditions at the boundaries y=±W that combine the discrete transparent boundary conditions (6.32) with source terms −ƒ_(j+1) ^((S))(±W) that are independent of the solutions f_(k)(y) for k≦j+1. Thus, $\begin{matrix} {{{{f_{j + 1}\left( {\pm W} \right)} \pm {\frac{1}{p_{0}}\frac{\partial}{\partial y}{f_{j + 1}\left( {\pm W} \right)}}} = {{{g_{j}\left( {\pm W} \right)}\quad{for}\quad j} = 0}},1,\ldots\quad,} & (7.26) \end{matrix}$  j=0,1, . . . ,  (7,26) where $\begin{matrix} {{g_{j}\left( {\pm W} \right)} = {{- {f_{j + 1}^{(S)}\left( {\pm W} \right)}} + {2p_{0}{\sum\limits_{k = 0}^{j}{\left( {- 1} \right)^{j - k}{{a_{k + 1}^{k}\left( {\pm W} \right)}.}}}}}} & (7.27) \end{matrix}$ The coefficients a_(k+1) ^(k) (±W) in this equation are computed recursively. First, a₁ ⁰(±W)=0. Second, for j=1,2, . . . , $\begin{matrix} {{{a_{k}^{j}\left( {\pm W} \right)} = {{{- \frac{1}{2}}{a_{k}^{j - 1}\left( {\pm W} \right)}} - {\frac{1}{2}{a_{k - 2}^{j - 1}\left( {\pm W} \right)}}}},{{for}\quad\underset{\underset{{{Skip}\quad{this}\quad{when}\quad j} = 1.}{︸}}{{k = 1},\ldots\quad,{j - 1}}},} & (7.28) \\ {{{a_{j}^{j}\left( {\pm W} \right)} = {{{\pm \quad\frac{1}{2}}{p_{0}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{f_{j - 1}\left( {\pm W} \right)}} + {\frac{\partial}{\partial y}{f_{j}\left( {\pm W} \right)}}} \right\rbrack}} - {a_{j}^{j - 1}\left( {\pm W} \right)} - {\frac{1}{2}{a_{j - 2}^{j - 1}\left( {\pm W} \right)}}}},} & (7.29) \\ {{{a_{j + 1}^{j}\left( {\pm W} \right)} = {{{\mp \frac{1}{2}}{p_{0}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{f_{j - 1}\left( {\pm W} \right)}} + {\frac{\partial}{\partial y}{f_{j}\left( {\pm W} \right)}}} \right\rbrack}} - {\frac{1}{2}{a_{j - 1}^{j - 1}\left( {\pm W} \right)}}}},} & (7.30) \end{matrix}$ where a⁻¹ ^(j)(±W)=0 for j=0,1, . . . ,  (7.31) a₀ ^(j)(±W)=0 for j=0,1, . . .  (7.32) If ƒ₀(y) is specified, then the N+1 equations in eqs.(7.14-7.19,7.26)) determine the polynomials ƒ_(j+1)(y) for j=0,1, . . . uniquely.

It is easy to evaluate ƒ_(j+1)(±W) from eq.(7.4) because T_(k)(±1)=(±1)^(k). $\begin{matrix} {{f_{j + 1}\left( {\pm W} \right)} = {\sum\limits_{k = 0}^{N}{\left( {\pm 1} \right)^{k}{f_{{j + 1},k}.}}}} & (7.33) \end{matrix}$ It is more convenient to evaluate σƒ_(j+1)(±W)/σy by differentiating eq.(7.4) directly than by using eq.(7.5) because T_(k)′(±1)=(±1)^(k+1)k². Thus, $\begin{matrix} {{{{\pm \frac{\partial}{\partial y}}f_{j + 1}} \pm (W)} = {\frac{1}{W}{\sum\limits_{k = 0}^{N}{\left( {\pm 1} \right)^{k}k^{2}{f_{{j + 1},k}.}}}}} & (7.34) \end{matrix}$ Now p₀W=e^(−iπ/4)2W(k₀Δx)^(1/2)=e^(−iπ/4)q^(1/2).  (7.35) Therefore, it follows from eqs.(7.26,7.33-7.35) that $\begin{matrix} {{{{\sum\limits_{k = 0}^{N}f_{{j + 1},k}} + {{\mathbb{e}}^{{\mathbb{i}}\quad{\pi/4}}q^{{- 1}/2}{\sum\limits_{k = 0}^{N}{k^{2}f_{{j + 1},k}}}}} = {g_{j}(W)}},} & (7.36) \\ {{{\sum\limits_{k = 0}^{N}{\left( {- 1} \right)^{k}f_{{j + 1},k}}} + {{\mathbb{e}}^{{\mathbb{i}}\quad{\pi/4}}q^{{- 1}/2}{\sum\limits_{k = 0}^{N}{\left( {- 1} \right)^{k}k^{2}f_{{j + 1},k}}}}} = {{g_{j}\left( {- W} \right)}.}} & (7.37) \end{matrix}$

It is easy to decouple the even-indexed Chebyshev coefficients from the odd-indexed Chebyshev coefficients in these equations. If we add eqs.(7.36,7.37), then we get an equation that involves only the even-indexed Chebyshev coefficients. $\begin{matrix} {{\sum\limits_{l = 0}^{N/2}{\left\lbrack {1 + {{\mathbb{e}}^{{\mathbb{i}}\quad{\pi/4}}{q^{{- 1}/2}\left( {2l} \right)}^{2}}} \right\rbrack f_{{j + 1},{2l}}}} = {{\frac{1}{2}\left\lbrack {{g_{j}(W)} + {g_{j}(W)}} \right\rbrack}.}} & (7.38) \end{matrix}$ If we subtract eq.(7.37) from eq.(7.36), then we get an equation that involves only the odd-indexed Chebyshev coefficients. $\begin{matrix} {{\sum\limits_{l = 1}^{N/2}{\left\lbrack {1 + {{\mathbb{e}}^{{\mathbb{i}}\quad{\pi/4}}{q^{{- 1}/2}\left( {{2l} - 1} \right)}^{2}}} \right\rbrack f_{{j + 1},{{2l} - 1}}}} = {{\frac{1}{2}\left\lbrack {{g_{j}(W)} - {g_{j}\left( {- W} \right)}} \right\rbrack}.}} & (7.39) \end{matrix}$

It is useful to express eqs.(7.14-7.19) and eqs.(7.38,7.39) as a pair of independent matrix equations of the form $\begin{matrix} {{\begin{bmatrix} A & u \\ v^{T} & \gamma \end{bmatrix}\begin{bmatrix} x \\ \alpha \end{bmatrix}} = {\begin{bmatrix} y \\ \beta \end{bmatrix}.}} & (7.40) \end{matrix}$ First, define $\begin{matrix} {\tau_{k} = \left\{ \begin{matrix} 0 & {{{{for}\quad k} = {N - 2}},{N - 3},} \\ {\frac{1}{4}\frac{1}{\left( {k + 1} \right)k}{iq}} & {{{{for}\quad k} = {N - 4}},\ldots\quad,2.} \end{matrix} \right.} & (7.41) \\ {s_{k} = \left\{ \begin{matrix} 1 & {{{{for}\quad k} = N},{N - 1},} \\ {1 - {\frac{1}{2}\frac{1}{k_{2} - 1}{iq}}} & {{{{for}\quad k} = {N - 2}},\ldots\quad,2.} \end{matrix} \right.} & (7.42) \\ {t_{k} = \left\{ \begin{matrix} {\frac{1}{4}\frac{1}{k\left( {k - 1} \right)}{iq}} & {{{{for}\quad k} = N},\ldots\quad,3,} \\ {\frac{1}{4}{iq}} & {{{for}\quad k} = 2.} \end{matrix} \right.} & (7.43) \end{matrix}$  b_(k)=1+e_(ix/4)q^(−1/2)k² for k=N, . . . ,0.  (7.44) The even-indexed Chebyshev coefficients ƒ_(j+1,2l) are determined by the matrix equation $\begin{matrix} {{{\begin{bmatrix} A_{0} & u_{0} \\ v_{0}^{T} & \gamma_{0} \end{bmatrix}\quad\begin{bmatrix} x_{0} \\ \alpha_{0} \end{bmatrix}} = \begin{bmatrix} y_{0} \\ \beta_{0} \end{bmatrix}},} & (7.45) \end{matrix}$ where $\begin{matrix} {{A_{0} = \begin{bmatrix} s_{N} & t_{N} & \quad & \quad & \quad & \quad & \quad \\ T_{N - 2} & s_{N - 2} & t_{N - 2} & \quad & \quad & \quad & \quad \\ \quad & {T_{N - 4}\quad} & s_{N - 4} & t_{N - 4} & \quad & \quad & \quad \\ \quad & \quad & ⋰ & ⋰ & ⋰ & \quad & \quad \\ \quad & \quad & \quad & ⋰ & ⋰ & ⋰ & \quad \\ \quad & \quad & \quad & \quad & r_{4} & s_{4} & t_{4} \\ \quad & \quad & \quad & \quad & \quad & r_{2} & s_{2} \end{bmatrix}},\quad{u_{0} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ \vdots \\ 0 \\ t_{2} \end{bmatrix}},} & (7.46) \\ {{v_{0} = \begin{bmatrix} b_{N} \\ b_{N - 2} \\ b_{N - 4} \\ \vdots \\ \vdots \\ b_{4} \\ b_{2} \end{bmatrix}},\quad{x_{0} = \begin{bmatrix} f_{{j + 1},N} \\ f_{{j + 1},{N - 2}} \\ f_{{j + 1},{N - 4}} \\ \vdots \\ \vdots \\ f_{{j + 1},4} \\ f_{{j + 1},2} \end{bmatrix}},\quad{y_{0} = \begin{bmatrix} h_{j,N} \\ h_{j,{N - 2}} \\ h_{j,{N - 4}} \\ \vdots \\ \vdots \\ h_{j,4} \\ h_{j,2} \end{bmatrix}},} & (7.47) \end{matrix}$ The rows of eq.(7.45) are eqs.(7.14,7.16), eq.(7.18) for k=2l (l=N/2−2, . . . ,2), eq.(7.19) and eq.(7.38).

The odd-indexed Chebyshev coefficients ƒ_(j+1,2l−1) are determined by the matrix equation $\begin{matrix} {{{\begin{bmatrix} A_{1} & u_{1} \\ v_{1}^{T} & \gamma_{1} \end{bmatrix}\quad\begin{bmatrix} x_{1} \\ \alpha_{1} \end{bmatrix}} = \begin{bmatrix} y_{1} \\ \beta_{1} \end{bmatrix}},} & (7.49) \end{matrix}$ where $\begin{matrix} {{A_{1} = \begin{bmatrix} s_{N - 1} & t_{N - 1} & \quad & \quad & \quad & \quad & \quad \\ \tau_{N - 3} & s_{N - 3} & t_{N - 3} & \quad & \quad & \quad & \quad \\ \quad & \tau_{N - 5} & s_{N - 5} & t_{N - 5} & \quad & \quad & \quad \\ \quad & \quad & ⋰ & ⋰ & ⋰ & \quad & \quad \\ \quad & \quad & \quad & ⋰ & ⋰ & ⋰ & \quad \\ \quad & \quad & \quad & \quad & \tau_{5} & s_{5} & t_{5} \\ \quad & \quad & \quad & \quad & \quad & \tau_{3} & s_{3} \end{bmatrix}},{u_{1} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ \vdots \\ 0 \\ t_{3} \end{bmatrix}},} & (7.50) \\ {{v_{1} = \begin{bmatrix} b_{N - 1} \\ b_{N - 3} \\ b_{N - 5} \\ \vdots \\ \vdots \\ b_{5} \\ b_{3} \end{bmatrix}},{x_{1}\begin{bmatrix} f_{{j + 1},{N - 1}} \\ f_{{j + 1},{N - 3}} \\ f_{{j + 1},{N - 5}} \\ \vdots \\ \vdots \\ f_{{j + 1},5} \\ f_{{j + 1},3} \end{bmatrix}},{y_{1} = \begin{bmatrix} h_{j,{N - 1}} \\ h_{j,{N - 3}} \\ h_{j,{N - 5}} \\ \vdots \\ \vdots \\ h_{j,5} \\ h_{j,3} \end{bmatrix}},} & (7.51) \end{matrix}$  γ₁=b₁, a₁=ƒ_(j+1,1), β₁=½|g_(j)(W)−g_(j)(−W)|.  (7.52) The rows of eq.(7.49) are eqs.(7.15,7.17), eq. (7.18) for k=2l−1(l=N/2−2, . . . ,2) and eq.(7.39).

We can solve eqs.(7.45,7.49) numerically by the Crout form of the block elimination method. The BEC algorithm to carry out this method has five steps. We outline them for eq.(7.40), which includes eqs.(7.45,7.49) as particular cases. Given y and β, we want to solve that equation for x and a.

(1) Solve the equation Aξ=u.  (7.53)

(2) Use the solution vector ξ to compute the divisor δ=γ−v^(T)·ξ.  (7.54)

(3) Solve the equation Aη=y.  (7.55)

(4) Use the solution vector η to compute the scalar a=(β−v^(T)·η)/δ.  (7.56)

(5) Compute the vector x=η−aξ.  (7.57) We can solve eqs.(7.53,7.55) with a method that uses an LU factorization of the matrix A. If A=A₀ or A=A₁, then the Crout LU factorization can be done numerically without pivoting and is stable.

Since s_(N)=1, we find that A₀=L₀U₀,  (7.58) where $\begin{matrix} {{L_{0} = \begin{bmatrix} 1 & \quad & \quad & \quad & \quad & \quad \\ \quad & \sigma_{N - 2} & \quad & \quad & \quad & \quad \\ \quad & \tau_{N - 4} & \sigma_{N - 4} & \quad & \quad & \quad \\ \quad & \quad & ⋰ & ⋰ & \quad & \quad \\ \quad & \quad & \quad & \tau_{4} & \sigma_{4} & \quad \\ \quad & \quad & \quad & \quad & \tau_{2} & \sigma_{2} \end{bmatrix}},{U_{0} = \begin{bmatrix} 1 & t_{N} & \quad & \quad & \quad & \quad \\ \quad & 1 & \tau_{N - 2} & \quad & \quad & \quad \\ \quad & \quad & 1 & \tau_{N - 4} & \quad & \quad \\ \quad & \quad & \quad & ⋰ & ⋰ & \quad \\ \quad & \quad & \quad & \quad & 1 & \tau_{4} \\ \quad & \quad & \quad & \quad & \quad & 1 \end{bmatrix}},} & (7.59) \end{matrix}$ and $\begin{matrix} \left. \begin{matrix} \begin{matrix} {{\sigma_{N - 2} = s_{N - 2}},} & {{\tau_{N - 2} = {t_{N - 2}/\sigma_{N - 2}}},} \end{matrix} \\ \begin{matrix} {{\sigma_{k} = {s_{k} - {\tau_{k}\tau_{k + 2}}}},} & {\tau_{k} = {t_{k}/\sigma_{k}}} & {{{{for}\quad k} = {2\quad{l\left( {{l = {{N/2} - 2}},\ldots\quad,2} \right)}}},} \end{matrix} \\ {\sigma_{2} = {s_{2} - {\tau_{2}{\tau_{4}.}}}} \end{matrix} \right\} & (7.60) \end{matrix}$ Since s_(N−1)=1, we find that A₁=L₁U₁,  (7.61) where $\begin{matrix} {{L_{1} = \begin{bmatrix} 1 & \quad & \quad & \quad & \quad & \quad \\ \quad & \sigma_{N - 3} & \quad & \quad & \quad & \quad \\ \quad & \tau_{N - 5} & \sigma_{N - 5} & \quad & \quad & \quad \\ \quad & \quad & ⋰ & ⋰ & \quad & \quad \\ \quad & \quad & \quad & \tau_{5} & \sigma_{5} & \quad \\ \quad & \quad & \quad & \quad & \tau_{3} & \sigma_{3} \end{bmatrix}},{U_{1} = \begin{bmatrix} 1 & t_{N - 1} & \quad & \quad & \quad & \quad \\ \quad & 1 & \tau_{N - 3} & \quad & \quad & \quad \\ \quad & \quad & 1 & \tau_{N - 5} & \quad & \quad \\ \quad & \quad & \quad & ⋰ & ⋰ & \quad \\ \quad & \quad & \quad & \quad & 1 & \tau_{5} \\ \quad & \quad & \quad & \quad & \quad & 1 \end{bmatrix}},} & (7.62) \end{matrix}$ and $\begin{matrix} \left. \begin{matrix} \begin{matrix} {{\sigma_{N - 3} = s_{N - 3}},} & {{\tau_{N - 3} = {t_{N - 3}/\sigma_{N - 3}}},} \end{matrix} \\ {{\sigma_{k} = {s_{k} - {\tau_{k}\tau_{k + 2}}}},{\tau_{k} = {{t_{k}/\sigma_{k}}\quad{for}}}} \\ {{k = {{2\quad l} - {1\left( {{l = {{N/2} - 2}},\ldots\quad,3} \right)}}},} \\ {\sigma_{3} = {s_{3} - {\tau_{3}{\tau_{5}.}}}} \end{matrix} \right\} & (7.63) \end{matrix}$

If A=A₀, then eqs.(7.53,7.55) are of the form A₀z₀=w₀. To solve this equation we introduce a vector p₀ and solve the first intermediate equation L₀p₀=w₀ for p₀ by forward substitution. Then we solve the second intermediate equation U₀z₀=p₀ for z₀ by backward substitution.

If A=A₁, then eqs.(7.53,7.55) are of the form A₁z₁=w₁. This equation is solved the same way. We introduce a vector p₁ and solve the intermediate equation L₁p₁=w₁ for p₁ by forward substitution. Then we solve the intermediate equation U₁z₁=p₁for z₁ by backward substitution.

Since the LU solution method for eqs.(7.53,7.55) is stable if A=A₀ or A=A₁, the BEC algorithm can be numerically unstable in finite-precision arithmetic only if the divisor δ is too small. However, if q=ο(N^(1/2)), then δ=ο(N^(1/2)) in each case. δ₀˜(2π)^(1/2)e^(−iπ/8)q^(1/4),  (7.64) δ₁˜(π/2)^(1/2)e^(−iπ/8)q^(1/4).  (7.65) Hence, if q=ο(N₂), then the BEC algorithm is stable. Each computer solution with the BEC algorithm should be followed by one round of iterative refinement. Solutions of eqs.(7.45,7.49) that are computed this way are accurate to machine precision. 8. Efficient evaluation of Chebyshev Interpolation polynomials.

Clenshaw's algorithm is a stable and accurate method to evaluate the Chebyshev interpolation polynomial $\begin{matrix} {{u(\xi)} = {\sum\limits_{k = 0}^{N}{u_{k}{T_{k}(\xi)}}}} & (8.1) \end{matrix}$ at a single argument ξ: u(ξ)=υ₀ where

(1) Start: υ_(N)=u_(N) and υ_(N−1)=u_(N−1)+2ξυ_(N).

(2) For k=N−2, . . . , 1: υ_(k)=u_(k)+2ξυ_(k+1)−υ_(k+2.)

(3) End: υ₀=u₀+ξυ₁−υ₂.

The arithmetic operation count for Clenshaw's algorithm is about 3N per evaluation. If we use this algorithm to evaluate u(ξ) at the N+1 points ξ_(n)=1−2n/N for n=0, . . . ,N,  (8.2) then the total operation count for these evaluations is ο(3N²). There are several alternative algorithms that also perform these evaluations but have total operation counts like ο(AN log N), where A is a constant. We shall describe an alternative algorithm that is based on Lagrange polynomial interpolation. If N is moderately large, say N>100, then it is more efficient to evaluate u(ξ) at all N+1 points ξ_(n) with this algorithm than with Clenshaw's algorithm.

If υ(θ)=u(cos θ), then N $\begin{matrix} {{v(\theta)} = {{\sum\limits_{l = 0}^{N}{u_{k}{\cos\left( {k\quad\theta} \right)}\quad{for}\quad 0}} \leq \theta \leq {\pi.}}} & (8.3) \end{matrix}$ Since ξ_(n)=cos φ_(n) where φ_(n)=2 arcsin{(n/N)^(1/2)} for n=0, . . . ,N,  (8.4) evaluating the polynomial u(ξ) at the regularly spaced points ξ_(n) is equivalent to evaluating υ(θ) at the irregularly spaced points φ_(n), which is what the alternative algorithm really does.

The first step of this algorithm is to compute the values of υ(θ) at the regularly spaced points θ_(j)=jπ/3N for j=0, . . . ,3N. A very efficient way to do these computations is to pad the spectral coefficients u_(n) in eq.(8.3) with 2N trailing zeros and then to take a cosine FFT of length 3N. The total operation count in this approach is ο(15N log N).

Next, let K be a positive integer. We shall discuss some choices for K after the algorithm has been described. The second step of the algorithm is to compute the forward differences of orders up to 2K that are associated with the points θ_(j). Let Δ⁰υ(θ_(j))=υ(θ_(j)) for j=0, . . . ,3N.  (8.5) For k=1, . . . ,2K the forward differences of order k are defined as follows. Δ^(k)υ(θ_(j))=Δ^(k−1)υ(θ_(j+1))−Δ^(k−1)υ(θ_(j)) for j=0, . . . ,3N−k.  (8.6) The operation count for computing all of these differences is ο(6KN). Note, however, that the subtractions in eq.(8.6) can be done in parallel for each k. Therefore, this step could be expedited on an array processor or on a computer that has a hardware implementation of vector arithmetic.

The third step of the algorithm is based on the Gaussian (forward) form of the Lagrange interpolation polynomials of even order. Let h=π/3N, and let Q_(2K)(θ_(j)+th) be the following polynomial in t of degree 2K. $\begin{matrix} {{Q_{2\quad K}\left( {\theta_{j} + {th}} \right)} = {{v\left( \theta_{j} \right)} + {t\quad\Delta\quad{v\left( \theta_{j} \right)}} + {\frac{1}{2}{t\left( {t - 1} \right)}\Delta^{2}{v\left( \theta_{j - 1} \right)}} + {\frac{1}{6}{t\left( {t^{2} - 1} \right)}\Delta^{3}{v\left( \theta_{j - 1} \right)}} + {\sum\limits_{k = 2}^{K - 1}{\prod\limits_{l = 1}^{k - 1}{\left( {t^{2} - l^{2}} \right)\left\lbrack {{\frac{t\left( {t - k} \right)}{\left( {2\quad k} \right)!}\Delta^{2\quad k}{v\left( \theta_{j - k} \right)}} + {\frac{t\left( {t^{2} - k^{2}} \right)}{\left( {{2\quad k} + 1} \right)!}\Delta^{{2\quad k} + 1}{v\left( \theta_{j - k} \right)}}} \right\rbrack}}} + {\prod\limits_{l = 1}^{K - 1}{\left( {t^{2} - l^{2}} \right)\frac{t\left( {t - K} \right)}{\left( {2\quad K} \right)!}\Delta^{2\quad K}{{v\left( \theta_{j - K} \right)}.}}}}} & (8.7) \end{matrix}$ Then υ(θ_(j)+th)=Q_(2K)(θ_(j)+th)+R_(2K)(θ_(j)+th),  (8.8) where the remainder R_(2K)(θ_(j)+th) can written in the form $\begin{matrix} {{R_{2\quad K}\left( {\theta_{j} + {th}} \right)} = {t{\prod\limits_{l = 1}^{K}{\left( {t^{2} - l^{2}} \right)\frac{h^{{2\quad K} + 1}}{\left( {{2\quad K} + 1} \right)!}\frac{\partial^{{2\quad K} + 1}}{\partial\theta^{{2\quad K} + 1}}{v(\vartheta)}}}}} & (8.9) \end{matrix}$ for some σ depending on t such that θ_(j−K)<σ<θ_(j+K). The polynomial Q_(2K)(θ_(j)+th) can be evaluated with about 4K operations as follows. $\begin{matrix} (1) & {{{Start}\text{:}\quad q_{K}} = {\frac{t - K}{2\quad K}\Delta^{2K}{{v\left( \theta_{j - K} \right)}.}}} \\ (2) & {{{{For}\quad k} = {K - 1}},\ldots\quad,{{1:q_{k}} = {{\frac{t - k}{2\quad k}\left\lbrack {{\Delta^{2\quad k}{v\left( \theta_{j - k} \right)}} + {\frac{t + k}{{2\quad k} + 1}\left( {{\Delta^{{2\quad k} + 1}{v\left( \theta_{j - k} \right)}} + q_{k + 1}} \right)}} \right\rbrack}.}}} \end{matrix}$

(3) End: Q_(2K)(θ_(j)+th)=υ(θ_(j))+t(Δυ(θ_(j))+q₁).

Now for n=0, . . . , N let j_(n) be the index j for which the absolute difference |φ_(n)−θ_(j)| is smallest, and let t_(n)=(φ_(n)−θ_(j) _(n) )/h. Clearly ${t_{n}} \leq {\frac{1}{2}.}$ Since φ₀=0=θ₀ and φ_(n)=π=θ_(3N), υ(φ₀)=υ(θ₀) and υ(φ_(n))=υ(θ_(3N)). If 0<n<N, then we use eq.(8.8) to approximate υ(φ_(n)) as follows. υ(φ_(n))=υ(θ_(j) _(n) +t_(n)h)≈Q_(2K)(θ_(j) _(n) +t_(n)h) for n=1, . . . ,N−1.  (8.10) The total operation count for computing these approximations is about 4KN. Note that in order for thereto be enough sample points θ_(j), the increment h=π/3N must be small enough that φ₁/h>K. A sufficient condition for there to be enough of these sample points is that N>(π²/36)K².

Let us estimate the interpolation errors R_(2K)(θ_(j) _(n) t_(n)h)=υ(φ_(n))−Q_(2K)(θ_(j) _(n) +t_(n)h)  (8.11) in the approximations (8.10). First, we consider the simple cases in which there is just one non-vanishing term in eq.(8.3). Let 0≦k≦N and assume that υ(θ)=cos(kθ). It can be shown that $\begin{matrix} {{{{R_{2K}\left( {\theta_{j_{n}} + {t_{n}h}} \right)}} \leq {\frac{1}{2}\left( \frac{\pi}{K} \right)^{\frac{1}{2}}{\left( {\frac{k}{N}\frac{\pi}{6}} \right)^{{2\quad K} + 1}\left\lbrack {1 + {O\left( K^{- 1} \right)}} \right\rbrack}\quad{if}}}\text{}{{v(\theta)} = {{\cos\left( {k\quad\theta} \right)}.}}} & (8.12) \end{matrix}$ For fixed k this bound is an exponentially decreasing function of increasing K because 0≦k/N≦1. For fixed K this bound is a rapidly decreasing function of decreasing k. Thus, the bound on the interpolation error is greatest if k=N. If K=25, however, eq.(8.12) implies double-precision accuracy even in this case because ½(π/25)^(1/2)(π/6)⁵¹=0.8×10⁻¹⁵.

In the general case we find that $\begin{matrix} {{{R_{2\quad K}\left( {\theta_{j_{n}} + {t_{n}h}} \right)}} \leq {\frac{1}{2}\left( \frac{\pi}{K} \right)^{\frac{1}{2}}{\sum\limits_{k = 1}^{N}{{u_{k}}{{\left( {\frac{k}{N}\frac{\pi}{6}} \right)^{{2\quad K} + 1}\left\lbrack {1 + {O\left( K^{- 1} \right)}} \right\rbrack}.}}}}} & (8.13) \end{matrix}$ In practice, the Chebyshev coefficients u_(k) eventually decrease exponentially in magnitude as the index k increases. Therefore, eq.(8.13) implies that the approximations in eq.(8.10) may have double-precision accuracy if the integer K is slightly less than 25, e.g., K≈20.

The total operation count for computing all the values u(ξ_(n))=υ(φ_(n)) with the alternative algorithm is ο(15N log N+10KN). This should be compared with the total operation count ο(3N²) for computing these values directly with Clenshaw's algorithm.

9. The decomposition method for transversely unbounded waveguides.

We describe a newly created numerical method for integrating a system of coupled parabolic equations that models one-way propagation of sound in a transversely unbounded waveguide. Referring now to FIG. 3, we assume that sound in the waveguide propagates mainly in the direction of increasing range 30, which is assigned the coordinate x. The waveguide is unbounded in the direction of cross range 32, which is assigned the coordinate y. This discussion is for slab geometry, but the development for spherical geometry is almost identical.

The horizontal variations in the properties of the fluid in this waveguide are confined between vertical planes at y=±W_(p), which are parallel to the direction of propagation. We assume that the fluid in the waveguide was originally in a horizontally stratified state but has been perturbed between these planes. The unperturbed fluid can have an internal horizontal interface that also has been perturbed between the planes. The perturbations are smooth and vanish at the planes. The waveguide remains horizontally stratified in the two semi-infinite sections where |y|>W_(p).

To make the transverse extent of our problem finite, we put transparent boundaries in the vertical planes at y=±W, where W>W_(p). These boundaries have no physical reality. Sound in the waveguide must penetrate the transparent boundaries without reflection because they are in the semi-infinite horizontally stratified sections of the waveguide.

For simplicity we assume that the reference wavenumbers k_(m) are positive constants. Thus, the acoustic pressure field is represented by modal amplitudes φ_(m) for m=1, . . . ,M that satisfy the following system of coupled parabolic equations for x>0. $\begin{matrix} {{{2\quad i\quad k_{m}\frac{\partial\varphi_{m}}{\partial x}} + \frac{\partial^{2}\varphi_{m}}{\partial y^{2}}} = {{\left( {k_{m}^{2} - \xi_{m}^{2}} \right)\varphi_{m}} + {\sum\limits_{n = 1}^{M}{\left\lbrack {{\left( {\alpha_{m\quad n} + {i\quad\beta_{m\quad n}k_{n}}} \right)\varphi_{n}} + {\beta_{m\quad n}\frac{\partial\varphi_{n}}{\partial x}} + {\gamma_{m\quad n}\frac{\partial\varphi_{n}}{\partial y}}} \right\rbrack e^{{- {i{({k_{m} - k_{n}})}}}x}}}}} & (9.1) \end{matrix}$ for m=1, . . . ,M. These equations are obtained from eq.(2.8) with σ_(m)=k_(m)x for each m.

Each modal wavenumber ξ_(m)(x, y) is constant in the half-planes 34 where |y|>W_(p), which correspond to the unperturbed horizontally stratified sections of the waveguide. We denote this constant by ξ_(m) ⁽⁰⁾. Note that it is the same for y<−W_(p) as for y>W_(p). We expect the reference wavenumber k_(m) to be close to ξ_(m) ⁽⁰⁾. In addition, the coupling coefficients a_(mn)(x,y),β_(mn)(x,y) and γ_(mn)(x,y) vanish for |y|>W_(p).

Sound enters the waveguide through an aperture 36 in the vertical plane at x=0. The aperture is confined between the lines at y=±W_(a) in this plane, where W_(a)<W_(p). The acoustic pressure field in this plane is represented by initial conditions on the modal amplitudes of the following form. $\begin{matrix} {{\varphi_{m}\left( {0,y} \right)} = \left\{ \begin{matrix} {S_{m}(y)} & {{{{if}\quad{y}} < W_{s}},} \\ 0 & {{{if}\quad{y}} > {W_{s}.}} \end{matrix} \right.} & (9.2) \end{matrix}$

Now let us formulate transparent boundary conditions for the modal amplitudes φ_(m) at the transparent boundaries 38 where y=±W. First, the coupled parabolic equations (9.1) decouple in the half-planes |y|>W_(p) as follows. For m=1, . . . ,M, $\begin{matrix} {{{2{ik}_{m}\frac{\partial\varphi_{m}}{\partial x}} + \frac{\partial^{2}\varphi_{m}}{\partial y^{2}}} = {{\left( {k_{m}^{2} - \left\{ \xi_{m}^{(0)} \right\}^{2}} \right)\varphi_{m}\quad{if}\quad{y}} > {W_{p}.}}} & (9.3) \end{matrix}$ We bring each of these equations into the form of a standard parabolic equation by defining the following phase-shifted modal amplitudes. {overscore (φ)}_(m)=φ_(m)e^(−1x) ^(m) ^(x) for m=1, . . . ,M,  (9.4) where each κ_(m) is a constant. If we let κ_(m)=({ξ_(m) ⁽⁰⁾}₂−k_(m) ²)/2 k_(m),  (9.5) then {overscore (φ)}_(m) satisfies the standard parabolic equation $\begin{matrix} {{{2{ik}_{m}\frac{\partial{\overset{\sim}{\varphi}}_{m}}{\partial x}} + \frac{\partial^{2}{\overset{\sim}{\varphi}}_{m}}{\partial y^{2}}} = {{0\quad{if}\quad{y}} > {W_{p}.}}} & (9.6) \end{matrix}$ Since W>W_(p), outgoing transparent boundary conditions can be imposed on the phase-shifted modal amplitudes at the transparent boundaries y=±W. For example, we get the following form of these conditions from eqs.(6.2,6.3). For m=1, . . . ,M, $\begin{matrix} {{{\overset{\sim}{\varphi}}_{m}\left( {x,{\pm W}} \right)} = {{{\mp {e^{i\quad{\pi/4}}\left( {2\quad\pi\quad k_{m}} \right)}^{{- 1}/2}}{\int_{0}^{x}{\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}\left( {u,{\pm W}} \right)}\left( {x - u} \right)^{{- 1}/2}{\mathbb{d}u}\quad{for}\quad x}}} > 0.}} & (9.7) \end{matrix}$

Let us reformulate eq.(9.1) in terms of the phase-shifted modal amplitudes. $\begin{matrix} {{{{2{ik}_{m}\frac{\partial{\overset{\sim}{\varphi}}_{m}}{\partial x}} + \frac{\partial^{2}{\overset{\sim}{\varphi}}_{m}}{\partial y^{2}}} = {{\left( {\left\{ \xi_{m}^{(0)} \right\}^{2} - \xi_{m}^{2}} \right){\overset{\sim}{\varphi}}_{m}} + {\sum\limits_{n = 1}^{M}{\left\lbrack {{\left( {a_{mn} + {i\quad\beta_{mn}{\overset{\sim}{k}}_{n}}} \right){\overset{\sim}{\varphi}}_{n}} + {\beta_{mn}\frac{\partial{\overset{\sim}{\varphi}}_{n}}{\partial x}} + {\gamma_{mn}\frac{\partial{\overset{\sim}{\varphi}}_{n}}{\partial y}}} \right\rbrack e^{{- {i{({{\overset{\sim}{k}}_{m} - {\overset{\sim}{k}}_{n}})}}}x}}}}},} & (9.8) \end{matrix}$ where the shifted reference wavenumbers {overscore (k)}_(m) are defined as {overscore (k)}_(m)=k_(m)+κ_(m)=ξ_(m) ⁽⁰⁾+(k_(m)−ξ_(m) ⁽⁰⁾)²/2k_(m) for m=1, . . . ,M.  (9.9) Note that {overscore (k)}_(m) is stationary at k_(m)=ξ_(m) ⁽⁰⁾, where it equals ξ_(m) ⁽⁰⁾ as well.

A sufficient condition that the system of coupled parabolic equations (9.1) have a unique solution for a given set of initial conditions (9.2) is that the reference matrix be positive definite and that the transparent boundary conditions (9.7) be satisfied.

Each {overscore (φ)}_(m) is the sum of a primary component {overscore (φ)}_(m) ⁽¹⁾ and a secondary component {overscore (φ)}_(m) ⁽²⁾ that are the solutions of the following two problems. I. Primary component. Let {overscore (φ)}_(m) ⁽¹⁾ be a solution of the inhomogeneous parabolic equation $\begin{matrix} {{{2{ik}_{m}\frac{\partial{\overset{\sim}{\varphi}}_{m}^{(1)}}{\partial x}} + \frac{\partial^{2}{\overset{\sim}{\varphi}}_{m}^{(1)}}{\partial y^{2}}} = {{\left( {\left\{ \xi_{m}^{(0)} \right\}^{2} - \xi_{m}^{2}} \right){\overset{\sim}{\varphi}}_{m}} + {\sum\limits_{n = 1}^{M}{\left\lbrack {{\left( {\alpha_{mn} + {i\quad\beta_{mn}{\overset{\sim}{k}}_{n}}} \right){\overset{\sim}{\varphi}}_{n}} + {\beta_{mn}\frac{\partial{\overset{\sim}{\varphi}}_{n}}{\partial x}} + {\gamma_{mn}\frac{\partial{\overset{\sim}{\varphi}}_{n}}{\partial y}}} \right\rbrack{\mathbb{e}}^{{- {{\mathbb{i}}{({{\overset{\sim}{k}}_{m} - {\overset{\sim}{k}}_{n}})}}}x}}}}} & (9.10) \end{matrix}$ for x>0 and |y|≦W that satisfies the same initial conditions that {overscore (φ)}_(m), satisfies, $\begin{matrix} {{{\overset{\sim}{\varphi}}_{m}^{(1)}\left( {0,y} \right)} = {{{\overset{\sim}{\varphi}}_{m}\left( {0,y} \right)} = \left\{ \begin{matrix} {S_{m}(y)} & {{{{if}\quad{y}} < W_{3}},} \\ 0 & {{{{if}\quad{y}} > W_{3}},} \end{matrix} \right.}} & (9.11) \end{matrix}$ and that is subject to the following Neumann boundary conditions. $\begin{matrix} {{\frac{\partial\quad}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(1)}\left( {x,{\pm W}} \right)}} = {{0\quad{for}\quad x} > 0.}} & (9.12) \end{matrix}$ Conditions (9.11,9.12) are compatible, and the solution {overscore (φ)}_(m) ⁽¹⁾ of this problem is unique. II. Secondary component. Let {overscore (k)}_(m) ⁽²⁾ be a solution of the standard parabolic equation $\begin{matrix} {{{2{ik}_{m}\frac{\partial{\overset{\sim}{\varphi}}_{m}^{(2)}}{\partial x}} + \frac{\partial^{2}{\overset{\sim}{\varphi}}_{m}^{(2)}}{\partial y^{2}}} = {{0\quad{for}\quad x} > {0\quad{and}\quad{y}} \leq W}} & (9.13) \end{matrix}$ that satisfies the homogeneous initial condition {overscore (φ)}_(m) ⁽²⁾(0,y)=0 for |y|≦W,  (9.14) and that is subject to the following inhomogeneous Neumann boundary conditions. $\begin{matrix} {{\frac{\partial\quad}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x,{\pm W}} \right)}} = {{\frac{\partial\quad}{\partial y}{{\overset{\sim}{\varphi}}_{m}\left( {x,{\pm W}} \right)}\quad{for}\quad x} > 0.}} & (9.15) \end{matrix}$ Conditions (9.14,9.15) are compatible, and the solution {overscore (φ)}_(m) ⁽²⁾ of this problem is unique.

It can be shown that {overscore (φ)}_(m)={overscore (φ)}_(m) ⁽¹⁾+{overscore (φ)}_(m) ⁽²⁾ for m=1, . . . ,M.  (9.16) This identity suggests the following basic decomposition method for finding the phase-shifted modal amplitudes {overscore (φ)}_(m). We regard each {overscore (φ)}_(m) as an unknown in problems I and II and solve these problems simultaneously for the triples {{overscore (φ)}_(m),{overscore (φ)}_(m) ⁽¹⁾,{overscore (φ)}_(m) ⁽²⁾)} subject to the constraints that {overscore (φ)}_(m)={overscore (φ)}_(m) ⁽¹⁾+{overscore (φ)}_(m) ⁽²⁾ and that {overscore (φ)}_(m) satisfy the transparent boundary conditions (9.7) for m=1, . . . ,M. This approach gives us the freedom to use different techniques for solving problems I and II numerically.

In practice, we solve discrete forms of problems I and II numerically. First, we develop a discrete form of problem I from an exact integral of eq.(9.10). Fix a range increment Δx and let x_(j)=jΔx for j=0,1, . . . . Define $\begin{matrix} {{{\overset{\sim}{\psi}}_{m}^{(1)} = {{k_{m}{\overset{\sim}{\psi}}_{m}^{(1)}} + {\left( {i/2} \right){\sum\limits_{n = 1}^{M}\quad{\theta_{mn}{\overset{\sim}{\psi}}_{n}{\mathbb{e}}^{{- {{\mathbb{i}}{({k_{m} - k_{n}})}}}x}}}}}},} & (9.17) \end{matrix}$ and $\begin{matrix} {{\overset{\_}{q}}_{m} = {{\left( {\left\{ \xi_{m}^{(0)} \right\}^{2} - \xi_{m}^{2}} \right){\overset{\sim}{\varphi}}_{m}} + {\sum\limits_{n = 1}^{M}{\left\lbrack {{\left( {\alpha_{mn} + {i{\overset{\_}{k}}_{m}\beta_{mn}} - \frac{\partial\beta_{mn}}{\partial x}} \right){\overset{\sim}{\varphi}}_{n}} + {\gamma_{mn}\frac{\partial{\overset{\sim}{\varphi}}_{n}}{\partial y}} + {\frac{i}{2k_{m}}\frac{\partial^{2}}{\partial y^{2}}\left( {\beta_{mn}{\overset{\sim}{\varphi}}_{n}} \right)}} \right\rbrack{e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x}.}}}}} & (9.18) \end{matrix}$ We can verify that eq.(9.10) is equivalent to the following equation. $\begin{matrix} {{{2{ik}_{m}\frac{\partial{\overset{\sim}{\psi}}_{m}^{(1)}}{\partial x}} + \frac{\partial^{2}{\overset{\sim}{\psi}}_{m}^{(1)}}{\partial y^{2}}} = {k_{m}{{\overset{\sim}{q}}_{m}.}}} & (9.19) \end{matrix}$ Note that {overscore (ψ)}_(m) ⁽¹⁾(x,y)=k_(m){overscore (φ)}_(m) ⁽¹⁾(x,y) and {overscore (q)}_(m)(x,y)=0 if |y|>W_(p)  (9.20) because ξ_(m)=ξ_(m) ⁽⁰⁾ and the coupling coefficients vanish for |y|>W_(p). Therefore, {overscore (ψ)}_(m) ⁽¹⁾(x, y) satisfies Neumann boundary conditions like eq.(9.12) at y=±W, and we can use the Fourier cosine transform ℑ to integrate eq.(9.19) from y=−W to y=W. Next, we use the integrating factor exp(iη²x/2k_(m)), where η is the spectral parameter, to integrate the transformed equation from x=x_(j) to x=x_(j+1). Applying ℑ⁻¹ to the result, we get ψ ~ m ( 1 ) ⁡ ( x j + 1 ) = ℱ - 1 ⁡ ( exp ⁢ { - i ⁢   ⁢ η 2 ⁢ Δ ⁢   ⁢ x / 2 / k m } ⁢ ℱ ⁢ ψ ~ m ( 1 ) ⁡ ( x j ) ) - ( i / 2 ) ⁢ ∫ x j x j + 1 ⁢ - 1 ⁢ ( - exp ⁢ { i ⁢   ⁢ η 2 ⁡ ( x j + 1 - x ) / 2 ⁢ k m } ⁢ ℱ ⁢   ⁢ q _ m ⁡ ( x ) ) ⁢ ⅆ x . ( 9.21 )

Now we introduce three numerical approximations in eq.(9.21) that lead to a system of linear equations for the values of {overscore (φ)}_(m) at discrete mesh points in the x, y plane. First, we use the trapezoidal rule to approximate the integral on the RHS of this equation. Dropping the ο((Δx)³) remainder in this approximation, we get {overscore (ψ)}_(m) ⁽¹⁾(x_(j+1))+(iΔx/4) {overscore (q)}_(m)(x_(j+1)) =ℑ⁻¹(exp{−iη²Δx/2k_(m)}ℑ[{overscore (ψ)}_(m) ⁽¹⁾(x_(j))−(iΔx/4{overscore (q)}_(m)(x_(j))]).  (9.22) This equation is inconsistent with eq.(9.19). However, it is consistent with a useful approximation to eq. (9.19), which is to regard {overscore (ψ)}_(m) ⁽¹⁾ as a solution of the standard parabolic equation $\begin{matrix} {{{2{ik}_{m}\frac{\partial{\overset{\sim}{\psi}}_{m}^{(1)}}{\partial x}} + \frac{\partial^{2}{\overset{\sim}{\psi}}_{m}^{(1)}}{\partial y^{2}}} = {{0\quad{for}\quad{y}} \leq W}} & (9.23) \end{matrix}$  for |y|≦W  (9.23) over the disjoint open intervals x_(j)<x<x_(j+1) for j=0,1, . . . that is subject to the following limit conditions at the range steps x_(j). {overscore (ψ)}_(m) ⁽¹⁾(x_(j)±0,y)={overscore (ψ)}_(m) ⁽¹⁾(x_(j),y)∓(iΔx/4){overscore (q)}_(m)(x_(j),y) for |y|<W.  (9.24) We can use these conditions to write eq.(9.22) as {overscore (ψ)}_(m) ⁽¹⁾(x_(j+1)−0)=ℑ⁻¹(exp{−iη²Δx/2k_(m)}ℑ[{overscore (ψ)}_(m) ⁽¹⁾(x_(j)+0)]),  (9.25) and we can combine them to get the following identities. $\begin{matrix} {{{{\overset{\sim}{\psi}}_{m}^{(1)}\left( {x_{j},y} \right)} = {\frac{1}{2}\left\lbrack {{{\overset{\sim}{\psi}}_{m}^{(1)}\left( {{x_{j} + 0},y} \right)} + {{\overset{\sim}{\psi}}_{m}^{(1)}\left( {{x_{j} - 0},y} \right)}} \right\rbrack}},} & (9.26) \\ {{{{\overset{\sim}{\psi}}_{m}^{(1)}\left( {{x_{j} + 0},y} \right)} - {{\overset{\sim}{\psi}}_{m}^{(1)}\left( {{x_{j} - 0},y} \right)}} = {\left( {{- i}\quad\Delta\quad{x/2}} \right){{{\overset{\_}{q}}_{m}\left( {x_{j},y} \right)}.}}} & (9.27) \end{matrix}$

Next, we fix a positive integer N, let Δy=2W/N and define the points y_(k)=kΔy−W for k=0,1, . . . ,N. Our second approximation is to replace the operators ℑ and ℑ⁻¹ in eq. (9.22) by the operators ℑ_(N) and ℑ_(N) ¹. Since these operators transform sequences to sequences, we must replace the functions {overscore (ψ)}_(m) ⁽¹⁾±(iΔx/4){overscore (q)}_(m) in this equation with the sequences of their values at the points y_(k). Let us write these sequences as {{overscore (ψ)}_(m) ⁽¹⁾(y_(k))±(iΔx/4){overscore (q)}_(m)(y_(k))}_(k=0) ^(N). Under these approximations, eq.(9.22) becomes {{overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),y_(k))+(iΔx/4){overscore (q)}_(m)(x_(j+1),y_(k))}_(k=0) ^(N)=ℑ⁻¹(exp{−iη²Δx/2k_(m)}ℑ_(N){overscore (ψ)}_(m) ⁽¹⁾(x_(j),y_(l))−(iΔx/4){overscore (q)}_(m)(x_(j),y_(l))}_(l=0) ^(N)  (9.28) The spectral parameter η in this equation is bounded because η=nπ/2W for n=0,1, . . . ,N. Hence, 0≦η≦π/Δy.

We must approximate the partial derivatives σ{overscore (φ)}_(n)/σy and σ²(β_(mn){overscore (φ)}_(n))/σ²y that appear in eq. (9.18) for {overscore (φ)}_(m) before we can evaluate the sequences in eq.(9.28). Since the variable y in this equation is restricted to the uniformly spaced points y_(k), our third approximation is to replace the partial derivatives at these points with central differences. At the lowest order of approximation, we use the three-point formulas $\begin{matrix} {{{\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{n}\left( y_{k} \right)}} = {{\frac{1}{2\quad\Delta\quad y}\left\lbrack {{{\overset{\sim}{\varphi}}_{n}\left( y_{k + 1} \right)} - {{\overset{\sim}{\varphi}}_{n}\left( y_{k - 1} \right)}} \right\rbrack} + {O\left( \left( {\Delta\quad y} \right)^{2} \right)}}},} & (9.29) \\ {{\frac{\partial^{2}}{\partial y^{2}}\left( {\beta_{mn}{\overset{\sim}{\varphi}}_{n}} \right)\left( y_{k} \right)} = {{\frac{1}{\left( {\Delta\quad y} \right)^{2}}\left\lbrack {{\left( {\beta_{mn}{\overset{\sim}{\varphi}}_{n}} \right)\left( y_{k + 1} \right)} - {2\left( {\beta_{mn}{\overset{\sim}{\varphi}}_{n}} \right)\left( y_{k} \right)} + {\left( {\beta_{mn}{\overset{\sim}{\varphi}}_{n}} \right)\left( y_{k - 1} \right)}} \right\rbrack} + {{O\left( \left( {\Delta\quad y} \right)^{2} \right)}.}}} & (9.30) \end{matrix}$

Now we combine these approximations with eq.(9.18). Rearranging terms in the result, we find that $\begin{matrix} {{{\overset{\_}{q}}_{m}\left( y_{k} \right)} = {{\left( {\left\{ \xi_{m}^{(0)} \right\}^{2} - {\xi_{m}^{2}\left( y_{k} \right)}} \right){{\overset{\sim}{\varphi}}_{m}\left( y_{k} \right)}} + {\sum\limits_{n = 1}^{M}{\left\lbrack {{\alpha_{mn}\left( y_{k} \right)} + {\left( {{i{\overset{\_}{k}}_{m}} - \frac{i}{{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right){\beta_{mn}\left( y_{k} \right)}} - {\frac{\partial}{\partial x}{\beta_{mn}\left( y_{k} \right)}}} \right\rbrack e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x}{{\overset{\sim}{\varphi}}_{n}\left( y_{k} \right)}}} + {\sum\limits_{n = 1}^{M}{\left( {{\frac{1}{2\Delta\quad y}{\gamma_{mn}\left( y_{k} \right)}} + {\frac{i}{2{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( y_{k + 1} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x}{{\overset{\sim}{\varphi}}_{n}\left( y_{k + 1} \right)}}} + {\sum\limits_{n = 1}^{M}{\left( {{{- \frac{1}{2\Delta\quad y}}{\gamma_{mn}\left( y_{k} \right)}} + {\frac{i}{2{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( y_{k - 1} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x}{{\overset{\sim}{\varphi}}_{n}\left( y_{k - 1} \right)}}} + {{O\left( \left( {\Delta\quad y} \right)^{2} \right)}.}}} & (9.31) \end{matrix}$ Therefore, neglecting terms of ο(Δx(Δy)²), we find that $\begin{matrix} {{{{\overset{\sim}{\psi}}_{m}^{(1)}\left( {x_{j + 1},y_{k}} \right)} + {\left( {i\quad\Delta\quad{x/4}} \right){{\overset{\_}{q}}_{m}\left( {x_{j + 1},y_{k}} \right)}}} = {{{\overset{\_}{k}}_{m}{{\overset{\sim}{\varphi}}_{m}^{(1)}\left( {x_{j + 1},y_{k}} \right)}} + {\left( {i\quad\Delta\quad{x/4}} \right)\left( {\left\{ \xi_{m}^{(0)} \right\}^{2} - {\xi_{m}^{2}\left( {x_{j + 1},y_{k}} \right)}} \right){{\overset{\sim}{\varphi}}_{m}\left( {x_{j + 1},y_{k}} \right)}} + {\left( {\frac{i}{2} - \frac{{\overset{\_}{k}}_{m}\Delta\quad x}{4} + \frac{\Delta\quad x}{4{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right){\sum\limits_{n = 1}^{M}{{\beta_{mn}\left( {x_{j + 1},y_{k}} \right)}e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}{{\overset{\sim}{\varphi}}_{n}\left( {x_{j + 1},y_{k}} \right)}}}} + {\left( {i\quad\Delta\quad{x/4}} \right){\sum\limits_{n = 1}^{M}{\left( {{\alpha_{mn}\left( {x_{j + 1},y_{k}} \right)} - {\frac{\partial}{\partial x}{\beta_{mn}\left( {x_{j + 1},y_{k}} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}{{\overset{\sim}{\varphi}}_{n}\left( {x_{j + 1},y_{k}} \right)}}}} + {\sum\limits_{n = 1}^{M}{\left( {{\frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}{\gamma_{mn}\left( {x_{j + 1},y_{k}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j + 1},y_{k + 1}} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}{{\overset{\sim}{\varphi}}_{n}\left( {x_{j + 1},y_{k + 1}} \right)}}} + {\sum\limits_{n = 1}^{M}{\left( {{{- \frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}}{\gamma_{mn}\left( {x_{j + 1},y_{k}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j + 1},y_{k - 1}} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}{{{\overset{\sim}{\varphi}}_{n}\left( {x_{j + 1},y_{k - 1}} \right)}.}}}}} & (9.32) \end{matrix}$

Next, we use the constraint that {overscore (φ)}_(m)={overscore (φ)}_(m) ⁽¹⁾+{overscore (φ)}_(m) ⁽²⁾ to replace the term k_(m){overscore (φ)}_(m) ⁽¹⁾(x₊₁,y_(k)) on the RHS of eq.(9.32) with k_(m){overscore (φ)}_(m)(x₊₁,y_(k)) by adding k_(m){overscore (φ)}_(m) ⁽²⁾(x₊₁,y_(k)) to each side of this equation. Thus, {overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),y_(k))+(iΔx/4){overscore (q)}_(m)(x_(j+1),y_(k))={overscore (ψ)}_(m)(x_(j+1)−0,y_(k))−k_(m){overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k))  (9.33) where $\begin{matrix} {{{\overset{\sim}{\Psi}}_{m}\left( {{x_{j + 1} - 0},y_{k}} \right)} = {{\left\lbrack {k_{m} + {\left( {i\quad\Delta\quad{x/4}} \right)\left( {\left\{ \xi_{m}^{(0)} \right\}^{2} - {\xi_{m}^{2}\left( {x_{j + 1},y_{k}} \right)}} \right)}} \right\rbrack{{\overset{\sim}{\varphi}}_{m}\left( {x_{j + 1},y_{k}} \right)}} + {\left( {\frac{i}{2} - \frac{{\overset{\_}{k}}_{m}\Delta\quad x}{4} + \frac{\Delta\quad x}{4{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right){\sum\limits_{n = 1}^{M}{\beta_{mn}\left( {x_{j + 1},y_{k}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}{{\overset{\sim}{\varphi}}_{n}\left( {x_{j + 1},y_{k}} \right)}}}} + {\left( {i\quad\Delta\quad{x/4}} \right){\sum\limits_{n = 1}^{M}{\left( {{\alpha_{mn}\left( {x_{j + 1},y_{k}} \right)} - {\frac{\partial}{\partial x}\beta_{mn}\left( {x_{j + 1},y_{k}} \right)}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}{{\overset{\sim}{\varphi}}_{n}\left( {x_{j + 1},y_{k}} \right)}}}} + {\sum\limits_{n = 1}^{M}{\left( {{\frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}{\gamma_{mn}\left( {x_{j + 1},y_{k}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j + 1},y_{k + 1}} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}{{\overset{\sim}{\varphi}}_{n}\left( {x_{j + 1},y_{k + 1}} \right)}}} + {\sum\limits_{n = 1}^{M}{\left( {{{- \frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}}{\gamma_{mn}\left( {x_{j + 1},y_{k}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j + 1},y_{k - 1}} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}{{{\overset{\sim}{\varphi}}_{n}\left( {x_{j + 1},y_{k - 1}} \right)}.}}}}} & (9.34) \end{matrix}$ Similarly, and to the same order of approximation, {overscore (ψ)}_(m) ⁽¹⁾(x_(j),y_(l))−(iΔx/4){overscore (q)}_(m)(x_(j),y_(l))={overscore (ψ)}_(m)(x_(j)+0,y_(l))−k_(m){overscore (ψ)}_(m) ⁽²⁾(x_(j),y_(l))  (9.35) where $\begin{matrix} {{{\overset{\sim}{\Psi}}_{m}\left( {{x_{j} - 0},y_{l}} \right)} = {{\left\lbrack {k_{m} - {\left( {i\quad\Delta\quad{x/4}} \right)\left( {\left\{ \xi_{m}^{(0)} \right\}^{2} - {\xi_{m}^{2}\left( {x_{j},y_{l}} \right)}} \right)}} \right\rbrack{{\overset{\sim}{\varphi}}_{m}\left( {x_{j},y_{l}} \right)}} + {\left( {\frac{i}{2} + \frac{{\overset{\_}{k}}_{m}\Delta\quad x}{4} - \frac{\Delta\quad x}{4{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right){\sum\limits_{n = 1}^{M}{\beta_{mn}\left( {x_{j},y_{l}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j}}{{\overset{\sim}{\varphi}}_{n}\left( {x_{j},y_{l}} \right)}}}} + {\left( {i\quad\Delta\quad{x/4}} \right){\sum\limits_{n = 1}^{M}{\left( {{\alpha_{mn}\left( {x_{j},y_{l}} \right)} - {\frac{\partial}{\partial x}{\beta_{mn}\left( {x,y_{l}} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}{{\overset{\sim}{\varphi}}_{n}\left( {x_{j},y_{l}} \right)}}}} + {\sum\limits_{n = 1}^{M}{\left( {{{- \frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}}{\gamma_{mn}\left( {x_{j},y_{l}} \right)}} + {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j},y_{l + 1}} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j}}{{\overset{\sim}{\varphi}}_{n}\left( {x_{j},y_{l + 1}} \right)}}} + {\sum\limits_{n = 1}^{M}{\left( {{\frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}{\gamma_{mn}\left( {x_{j},y_{l}} \right)}} - {\frac{\Delta\quad x}{8{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{mn}\left( {x_{j},y_{l - 1}} \right)}}} \right)e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j}}{{{\overset{\sim}{\varphi}}_{n}\left( {x_{j},y_{l - 1}} \right)}.}}}}} & (9.36) \end{matrix}$

We shall assume that Δy<W−W_(p). Under this assumption the coupling coefficients a_(mn)(x,y_(k)), β_(mn)(x,y_(k)) and γ_(mn)(x,y_(k)) vanish for x>0 and k=1, N−1. This eliminates the boundary values {overscore (σ)}_(m)(x_(j+1),±W) from the RHS of eq.(9.34), which makes it possible for us to find the internal values {overscore (σ)}_(m)(x_(j+1),y_(k)) (k=1, . . . ,N−1) from this equation. Let us combine eqs.(9.28,9.33,9.35) as follows. {{overscore (Ψ)}_(m)(x_(j+1)−0,y_(k))−k_(m){overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k))}_(k=0) ^(N) =ℑ_(N) ⁻¹(exp{−iη²Δx/2k_(m)}ℑ_(N){{overscore (ψ)}_(m)(x_(j)+0,y_(l))−k_(m){overscore (ψ)}_(m) ⁽²⁾(x_(j),y_(l))}_(l=0) ^(N)).  (9.37) If the values {overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k))(k=1, . . . ,N−1) at the range step x_(j+1) are known, then eqs.(9.34,9.37) determine a system of linear equations for the values {overscore (ψ)}_(m)(x_(j+1),y_(k))(k=1, . . . ,N−1) at the range step x_(j+1). Let us put these equations in matrix form. First, define the column vectors {overscore (ψ)}_(k)=[{overscore (ψ)}₁(x_(j+1),y_(k)), . . . ,{overscore (ψ)}_(M)(x_(j+1),y_(k))]^(T)  (9.38) {overscore (Ψ)}_(k) ⁻=[{overscore (Ψ)}₁(x_(j+1)−0,y_(k)), . . . ,{overscore (Ψ)}_(M)(x_(j+1)−0,y_(k)]^(T),  (9.39) for k=1, . . . ,N−1. Next, define the M×M matrices A_(kk) and A_(kk+1) as follows. $\begin{matrix} {\left( A_{k\quad k} \right)_{m\quad n} = \left\{ \begin{matrix} {{k_{n} + {\left( {i\quad\Delta\quad{x/4}} \right)\left( {\left\{ \xi_{n}^{(0)} \right\}^{2} - {\xi_{n}^{2}\left( {x_{j + 1},y_{k}} \right)} + {\alpha_{n\quad n}\left( {x_{j + 1},y_{k}} \right)}} \right)}},} & \left( {{{if}\quad m} = n} \right) \\ \left\lbrack {{\left( {\frac{i}{2} - \frac{{\overset{\_}{k}}_{m}\Delta\quad x}{4} + \frac{\Delta\quad x}{4\quad{k_{m}\left( {\Delta\quad y} \right)}^{2}}} \right){\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}} +} \right. & \quad \\ {{\left. {\left( {i\quad\Delta\quad{x/4}} \right)\left( {{\alpha_{m\quad n}\left( {x_{j + 1},y_{k}} \right)} - {\frac{\partial}{\partial x}{\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}}} \right)} \right\rbrack e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}},} & \left( {{{if}\quad m} \neq n} \right) \end{matrix} \right.} & (9.40) \\ {\left( A_{{k\quad k} \pm 1} \right)_{m\quad n} = {\left( {{{\pm \frac{i\quad\Delta\quad x}{8\Delta\quad y}}{\gamma_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}} - {\frac{\Delta\quad x}{8\quad{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{m\quad n}\left( {x_{j + 1},y_{k \pm 1}} \right)}}} \right){e^{{- {i{({{\overset{\_}{k}}_{m} - {\overset{\_}{k}}_{n}})}}}x_{j + 1}}.}}} & (9.41) \end{matrix}$

The assumption that Δy<W−W_(p) eliminates the boundary values {overscore (ψ)}_(m)(x_(j+1), ±W) from the RHS of eq.(9.34). Thus, it follows from this equation that $\begin{matrix} {{\begin{bmatrix} A_{11} & A_{12} & \quad & \quad & \quad \\ A_{21} & A_{22} & A_{23} & \quad & \quad \\ \quad & ⋰ & ⋰ & ⋰ & \quad \\ \quad & \quad & A_{N - {2N} - 3} & A_{N - {2N} - 2} & A_{N - {2N} - 1} \\ \quad & \quad & \quad & A_{N - {1N} - 2} & A_{N - {1N} - 1} \end{bmatrix} \cdot \begin{bmatrix} {\overset{\sim}{\varphi}}_{1} \\ {\overset{\sim}{\varphi}}_{2} \\ \vdots \\ {\overset{\sim}{\varphi}}_{N - 2} \\ {\overset{\sim}{\varphi}}_{N - 1} \end{bmatrix}} = {\begin{bmatrix} {\overset{\sim}{\Psi}}_{\overset{\_}{1}} \\ {\overset{\sim}{\Psi}}_{\overset{\_}{2}} \\ \vdots \\ {\overset{\sim}{\Psi}}_{\overset{\_}{N} - 2} \\ {\overset{\sim}{\Psi}}_{\overset{\_}{N} - 1} \end{bmatrix}.}} & (9.42) \end{matrix}$ If we use eq.(9.37) and the values {overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k))(k=1, . . . ,N−1) to evaluate the vectors {overscore (Ψ)}_(k) ⁻ on the RHS of this equation, then we can solve it for the internal solution vectors {overscore (ψ)}_(k).

We can solve eq.(9.42) numerically by Gaussian elimination with partial pivoting (GEPP) followed by one round of iterative refinement. The block tri-diagonal system matrix on the LHS of eq.(9.42) is also a band matrix of order M(N−1) that has 2M−1 lower diagonals and 2M−1 upper diagonals. Since M<<N in practice, it is most efficient to do the GEPP computations with subroutines that take advantage of band storage.

The system matrix on the LHS of eq. (9.42) is nonsingular if it is strictly diagonally dominant by columns. We see from eqs.(9.40,9.41) that this condition requires the following strict inequalities to be satisfied for the block indices k=1, . . . ,N−1:

|k_(n)+(iΔx/4)({ξ_(n) ⁽⁰⁾}²−ξ_(n) ²(x_(j+1)y_(k))+a_(nn)(x_(j+1),y_(k)))| $\begin{matrix} \begin{matrix} {{{k_{n} + {\left( {i\quad\Delta\quad{x/4}} \right)\left( {\left\{ \xi_{n}^{(0)} \right\}^{2} - {\xi_{n}^{2}\left( {x_{j + 1},y_{k}} \right)} + {\alpha_{n\quad n}\left( {x_{j + 1},y_{k}} \right)}} \right)}}} >} \\ {\sum\limits_{\underset{({m \neq n})}{m = 1}}^{M}{{{\left( {\frac{i}{2} - \frac{{\overset{\_}{k}}_{m}\Delta\quad x}{4} + \frac{\Delta\quad x}{4\quad{k_{m}({\Delta y})}^{2}}} \right){\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}} +}}} \\ {{{\left( {i\quad\Delta\quad{x/4}} \right)\left( {{\alpha_{m\quad n}\left( {x_{j + 1},y_{k}} \right)} - {\frac{\partial}{\partial x}{\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}}} \right)}} +} \\ {{\sum\limits_{m = 1}^{M}{{{\frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}{\gamma_{m\quad n}\left( {x_{j + 1},y_{k - 1}} \right)}} - {\frac{\Delta\quad x}{8\quad{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}}}}} +} \\ {\sum\limits_{m = 1}^{M}{{{{- \frac{i\quad\Delta\quad x}{8\quad\Delta\quad y}}{\gamma_{m\quad n}\left( {x_{j + 1},y_{k + 1}} \right)}} - {\frac{\Delta\quad x}{8\quad{k_{m}\left( {\Delta\quad y} \right)}^{2}}{\beta_{m\quad n}\left( {x_{j + 1},y_{k}} \right)}}}}} \\ {{{{for}\quad n} = 1},\ldots\quad,{M.}} \end{matrix} & (9.43) \end{matrix}$  for n=1, . . . ,M.  (9.43)

If this condition is satisfied, then Gaussian elimination requires no pivoting and is stable.

Note the factors {overscore (k)}_(m)Δx, Δx/Δy and Δx/k_(m)(Δy)² that appear in this equation. We recommend choosing the increments Δx and Δy so that k_(m)Δx≦ο(1) and k_(m)Δy≧ο(1). These restrictions imply that the ratios Δx/Δy and Δx/k_(m)(Δy)² are at most ο(1).

If Δy<W−W_(p), then we can find the boundary values {overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),±W) of the primary component from eq.(9.37) without the values {overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k)). In this case, it follows from eq.(9.32) that {overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),±W)+(iΔx/4){overscore (q)}_(m)(x_(j+1),±W)=k_(m){overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),±W). Thus, eqs.(9.33,9.37) imply that k m ⁢ ψ ~ m ( 1 ) ⁡ ( x j + 1 , ± W ) = N - 1 ( vh ± w ) ⁢ ( exp ⁢ { - i ⁢   ⁢ η 2 ⁢ Δ ⁢   ⁢ x / 2 ⁢ k m } ⁢ N ⁢ { Ψ ~ m ⁡ ( x j + 0 , yl ) - k m ⁢ ψ ~ m ( 2 ) ⁡ ( x j , yn ) } t = 0 N ) . ( 9.44 ) Note that it is not necessary to compute the whole inverse DFT to evaluate the RHS of this equation because we need only the values from the inverse DFT at the endpoints y_(k)=±W.

This completes our development of the discrete form of problem I and the method for solving it. Next, we develop a discrete form of problem II. We begin by integrating eq.(9.13) from x=x_(j) to x=x_(j+1) with the trapezoidal rule. Dropping the ο((Δx)³) remainder in this approximation, we get the semi-discrete parabolic equation $\begin{matrix} {{{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j + 1},y} \right)} - {\overset{\sim}{\varphi_{m}^{(2)}}\left( {x_{j},y} \right)}} = {{{\left( {i\quad\Delta\quad{x/4}\quad k_{m}} \right)\left\lbrack {{\frac{\partial^{2}}{\partial y^{2}}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j},y} \right)}} + {\frac{\partial^{2}}{\partial y^{2}}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j + 1},y} \right)}}} \right\rbrack}\quad{for}\quad{y}} \leq {W.}}} & (9.45) \end{matrix}$ We adopt conditions (9.14,9.15) without change. In particular, we require that the inhomogeneous Neumann boundary conditions (9.15) be satisfied at each range step. $\begin{matrix} {{\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j + 1},{\pm W}} \right)}} = {\frac{\partial}{\partial y}{{{\overset{\sim}{\varphi}}_{m}\left( {x_{j + 1},{\pm W}} \right)}.}}} & (9.46) \end{matrix}$

In order to formulate transparent boundary conditions for the discretized amplitude {overscore (ψ)}_(m), we assume that it satisfies the semi-discrete parabolic equation that corresponds to eq.(9.6). $\begin{matrix} {{{{\overset{\sim}{\varphi}}_{m}\left( {x_{j + 1},y} \right)} - {{\overset{\sim}{\varphi}}_{m}\left( {x_{j},y} \right)}} = {{{\left( {i\quad\Delta\quad{x/4}\quad k_{m}} \right)\left\lbrack {{\frac{\partial}{\partial y^{2}}{{\overset{\sim}{\varphi}}_{m}\left( {x_{j},y} \right)}} + {\frac{\partial^{2}}{\partial y^{2}}{{\overset{\sim}{\varphi}}_{m}\left( {x_{j + 1},y} \right)}}} \right\rbrack}\quad{for}\quad{y}} > {W_{p}.}}} & (9.47) \end{matrix}$ Under this assumption, the discretized amplitude {overscore (ψ)}_(m) satisfies discrete transparent boundary conditions like eq.(6.31). We can adapt the algorithm described by eqs.(6.26-6.30) to formulate these conditions as follows. $\begin{matrix} {{{{{\overset{\sim}{\varphi}}_{m}\left( {x_{j + 1},{\pm W}} \right)} \pm {\frac{1}{p_{m}}\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}\left( {x_{j + 1},{\pm W}} \right)}}} = {{2\quad p_{m}{\sum\limits_{k = 0}^{j}{\left( {- 1} \right)^{j - k}{a_{{k + 1},m}^{k}\left( {\pm W} \right)}\quad{for}\quad j}}} = 0}},{1\ldots}\quad,} & (9.48) \end{matrix}$ where p_(m)=e^(−iπ/4)2(k_(m)/Δx)^(1/2)  (9.49) and the coefficients a_(k+1,m) ^(k)(±W) are computed recursively. First, a_(1,m) ⁰(±W)=0. Second, for j=1,2, . . . , $\begin{matrix} {{{a_{k,m}^{j}\left( {\pm W} \right)} = {{{- \frac{1}{2}}{a_{k,m}^{j - 1}\left( {\pm W} \right)}} - {\frac{1}{2}{a_{{k - 2},m}^{j - 1}\left( {\pm W} \right)}}}},{{for}\quad\underset{\underset{{{Skip}\quad{this}\quad{when}\quad j} = 1.}{︸}}{\left( {{k = 1},\ldots\quad,{j - 1}} \right)}},} & (9.50) \\ {{{a_{j,m}^{j}\left( {\pm W} \right)} = {{{\pm \frac{1}{2}}{p_{m}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}\left( {x_{j - 1},{\pm W}} \right)}} + {\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}\left( {x_{j},{\pm W}} \right)}}} \right\rbrack}} - {a_{j,m}^{j - 1}\left( {\pm W} \right)} - {\frac{1}{2}{a_{{j - 2},m}^{j - 1}\left( {\pm W} \right)}}}},} & (9.51) \\ {{{a_{{j + 1},m}^{j}\left( {\pm W} \right)} = {{{\mp \frac{1}{2}}{p_{m}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}\left( {x_{j - 1},{\pm W}} \right)}} + {\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}\left( {x_{j},{\pm W}} \right)}}} \right\rbrack}} - {\frac{1}{2}{a_{{j - 1},m}^{j - 1}\left( {\pm W} \right)}}}},} & (9.52) \end{matrix}$ where a_(−1,m) ^(j)(±W)=0, a_(0,m) ^(j)(±W)=0 for j=0,1, . . .  (9.53)

Now we combine the constraints {overscore (ψ)}_(m)(x_(j+1),±W)={overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),±W)+{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),±W) and eq.(9.49) with eq.(9.48) to develop new boundary conditions for {overscore (ψ)}_(m) ⁽²⁾. Note that the boundary values of the primary component {overscore (ψ)}_(m) ⁽¹⁾ appear as source terms in these conditions. $\begin{matrix} {{{{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j + 1},{\pm W}} \right)} \pm {\frac{1}{p_{m}}\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j + 1},{\pm W}} \right)}}} = {{{- {{\overset{\sim}{\varphi}}_{m}^{(1)}\left( {x_{j + 1},{\pm W}} \right)}} + {2\quad p_{m}{\sum\limits_{k = 0}^{j}{\left( {- 1} \right)^{j - k}{a_{{k + 1},m}^{k}\left( {\pm W} \right)}\quad{for}\quad j}}}} = 0}},1,{\ldots\quad.}} & (9.54) \end{matrix}$ We also combine condition (9.46) with recursion formulas (9.51,9.52) to reformulate them in terms of {overscore (104 )}_(m) ⁽²⁾. $\begin{matrix} {{{a_{j,m}^{j}\left( {\pm W} \right)} = {{{\pm \frac{1}{2}}{p_{m}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j - 1},{\pm W}} \right)}} + {\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j},{\pm W}} \right)}}} \right\rbrack}} - {a_{j,m}^{j - 1}\left( {\pm W} \right)} - {\frac{1}{2}{a_{{j - 2},m}^{j - 1}\left( {\pm W} \right)}}}},} & (9.55) \\ {{a_{{j + 1},m}^{j}\left( {\pm W} \right)} = {{{\mp \frac{1}{2}}{p_{m}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{{\hat{\varphi}}_{m}^{(2)}\left( {x_{j - 1},{\pm W}} \right)}} + {\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j},{\pm W}} \right)}}} \right\rbrack}} - {\frac{1}{2}{{a_{{j - 1},m}^{j - 1}\left( {\pm W} \right)}.}}}} & (9.56) \end{matrix}$ Therefore, in principle the discrete form of problem II consists of

-   -   (1) the semi-discrete parabolic equation (9.45), subject to     -   (2) the initial condition (9.14), and     -   (3) the boundary conditions (9.54), with     -   (4) the recursion formulas (9.50,9.53,9.55,9.56).         However, in practice we would use the fixed-length truncated         algorithm described by eqs.(6.72-6.76) to evaluate the boundary         conditions for this problem. Let L_(J) denote the length of the         algorithm, which could be chosen according to criterion. (6.82)         for calculations in double-precision arithmetic. First, we use         eq.(6.78) to modify the boundary conditions (9.54) as follows.         $\begin{matrix}         {{{{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j + 1},{\pm W}} \right)} \pm {\frac{1}{p_{m}}\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j + 1},{\pm W}} \right)}}} = {{{- {{\overset{\sim}{\varphi}}_{m}^{(1)}\left( {x_{j + 1},{\pm W}} \right)}} + {2\quad p_{m}{\sum\limits_{k = 0}^{j}{\left( {- 1} \right)^{j - k}{\alpha_{L_{J},m}^{k}\left( {\pm W} \right)}\quad{for}\quad j}}}} = 0}},1,{\ldots\quad.}} & (9.57)         \end{matrix}$         Second, we use eqs.(6.72-6.76) to replace         eqs.(9.50,9.53,9.55,9.56) with the following recursion equations         for the coefficients a_(L) _(J) _(m) ^(k)(±W). $\begin{matrix}         {{{\alpha_{l,m}^{j}\left( {\pm W} \right)} = {{{- \frac{1}{2}}{\alpha_{{l + 1},m}^{j - 1}\left( {\pm W} \right)}} - {\frac{1}{2}{\alpha_{{l - 1},m}^{j - 1}\left( {\pm W} \right)}}}}{{{{for}\quad l} = 1},\ldots\quad,{L_{J} - 2},}} & (9.58) \\         {{{\alpha_{{L_{J} - 1},m}^{j}\left( {\pm W} \right)} = {{{\pm \frac{1}{2}}{p_{m}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j - 1},{\pm W}} \right)}} + {\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j},{\pm W}} \right)}}} \right\rbrack}} - {\alpha_{L_{J},m}^{- 1}\left( {\pm W} \right)} - {\frac{1}{2}{\alpha_{{L_{J} - 2},m}^{j - 1}\left( {\pm W} \right)}}}},} & (9.59) \\         {{{\alpha_{L_{J},m}^{j}\left( {\pm W} \right)} = {{{\mp \frac{1}{2}}{p_{m}^{- 2}\left\lbrack {{\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j - 1},{\pm W}} \right)}} + {\frac{\partial}{\partial y}{{\overset{\sim}{\varphi}}_{m}^{(2)}\left( {x_{j},{\pm W}} \right)}}} \right\rbrack}} - {\frac{1}{2}{\alpha_{{L_{J} - 1},m}^{j - 1}\left( {\pm W} \right)}}}},} & (9.60)         \end{matrix}$         where         a_(l,m) ⁰(±W)=0 for l=1, . . . ,L_(J), a_(0,m) ^(j)(±W)=0 for         j=0,1, . . .  (9.61)

This completes the development of the discrete form of problem II. A sufficient condition for it to be stable is that $\begin{matrix} {{\frac{\Delta\quad x}{{k_{0}\left( {\Delta\quad y} \right)}^{2}} < \frac{2}{\pi}},} & (9.62) \end{matrix}$ where k₀=min k_(m)(m=1, . . . ,M). Numerical integration of eq.(9.45) can introduce numerical artifacts in the solutions of the discrete form of problem II. However, these artifacts will be suppressed if the spatial frequency limit of the integration method does not exceed the following bound. λ_(c)=3^(−1/4)2(k₀/Δx)^(1/2).  (9.63) This criterion imposes a simple constraint on the N-th order Chebyshev tau integration method. Since the spatial frequency limit of the N-th order Chebyshev tau integration method is 2/Δy, this limit is less than λ_(c) if $\begin{matrix} {\frac{\Delta\quad x}{{k_{0}\left( {\Delta\quad y} \right)}^{2}} < {3^{{- 1}/2}.}} & (9.64) \end{matrix}$ This condition is only slightly more restrictive than condition (9.62).

Referring now to FIG. 4, we outline the operation of the preferred embodiment. We begin the process in block 40 by computing the modal wavenumbers and the coupling coefficients at the mesh points (x_(j),y_(k))(k=0,1, . . . ,N) from environmental data. They can be computed for all j=0,1, . . . and stored for subsequent use or computed as needed, beginning at j=0. In block 42 we select the initial values {overscore (ψ)}_(m)(x₀,y_(k))(j=0 and k=0,1, . . . ,N) for the source field in the waveguide at the plane x=x₀. Now we go to block 44. This block is the start of the general cycle of steps for solving problems I and II simultaneously. We must compute the values {overscore (ψ)}_(m)(x_(j+1),y_(k)) and {overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k)) for m=1, . . . ,M and k=0,1, . . . ,N at an arbitrary range step x_(j)+1 from the values {overscore (ψ)}_(m)(x_(j),y_(k)) and {overscore (ψ)}_(m) ⁽²⁾(x_(j),y_(l)) for m=1, . . . ,M and l=0,1, . . . ,N at the preceding range step x_(j).

(1) In block 46, for each m we compute the boundary values {overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),±W) of the primary component from eq.(9.44). These values are used in the next step.

(2) In block 48, for each m we integrate the semi-discrete parabolic equation (9.45) subject to the boundary conditions (9.57), which we evaluate with eqs.(9.58-9.61). We use the Chebyshev tau method to do the integration numerically. This method transforms the integration to an algebraic problem for the Chebyshev coefficients of {overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y) that we solve with the BEC algorithm. It is then necessary to compute the values {overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k)) for k=0,1, . . . ,N from these coefficients. Depending on the size of N, we use either Clenshaw's algorithm or an alternative algorithm to compute them. We use the internal values {overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k))(k=1, . . . ,N−1) in the next step. We compute the differential boundary values σ{overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),±W)/σy and save them for use by eqs.(9.59,9.60) at the following two range steps, x_(j+2) and x_(j+3). Having found {overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),±W) in the previous step and {overscore (ψ)}_(m) ⁽²⁾(x_(j+1),±W) in the present step, we compute the boundary values {overscore (ψ)}_(m)(x_(j+1),±W) from the constraints that {overscore (ψ)}_(m)(x_(j+1),±W)={overscore (ψ)}_(m) ⁽¹⁾(x_(j+1),±W)+{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),±W).

(3) In block 50, we solve eq.(9.42) numerically by GEPP for the internal values {overscore (ψ)}_(m)(x_(j+1),y_(k)) (k=1, . . . ,N−1). This completes the general cycle of steps for solving problems I and II.

In block 52 we use the modal amplitudes to evaluate the eigenfunction expansion of the complex pressure at desired depths in the waveguide at range x_(j+1) and cross ranges y_(k)(k=0,1, . . . ,N). We predict the physical acoustic pressure from these computations. Finally, we increase j to j+1 in block 54 and return to block 44 for the next cycle of computations. This completes the outline.

The discretized solutions {overscore (ψ)}_(m) turn out to be inaccurate at very large range from the source due to small inaccuracies in the discretized solutions {overscore (ψ)}_(m) ⁽²⁾. However, an extrapolation method can be applied to the discretized solutions {overscore (ψ)}_(m) ⁽²⁾ to compensate for the inaccuracies in them under conditions of practical interest. Before describing this method, let us modify the notation for these solutions to indicate their dependence on the increment Δx. {overscore (ψ)}_(m) ⁽²⁾(x_(j),y)→{overscore (ψ)}_(m) ⁽²⁾(x_(j),y_(i)Δx)  (9.65) A simple way to improve the accuracy of these approximations is to reduce the step size Δx. Since they are even functions of Δx, it is efficient to combine this approach with Richardson extrapolation.

We shall discuss a procedure that uses Richardson extrapolation to estimate {overscore (ψ)}_(m) ⁽²⁾(x_(j),y) from the computed values of {overscore (ψ)}_(m) ⁽²⁾(x_(j),y;Δx/υ) for υ=1,2,3,4 if Δx=ο(π/k₀). First, we introduce some new notation. Define the interpolated range steps x_(μ) ^((υ)) as follows. x_(μ) ^((υ))=μΔx/y for μ=1,2, . . . and υ=1,2, . . .  (9.66) Note that x_(υj) ^((υ)=x) _(j) for all j and υ.

The procedure requires that the discrete form of problem II be solved four times over each interval x_(j)≦x≦x_(j+1) using fractional steps Δx/υ for υ=1,2,3,4. For each υ we must replace Δx in eqs.(9.45,9.49) with Δx/υ and solve this problem at the interpolated range steps x_(υj+μ) ^((υ)) for μ=1, . . . ,υ. If we use the Chebyshev tau method to do the cross range integration in eq.(9.45), then we can do these calculations in terms of the Chebyshev expansion coefficients of the solutions {overscore (ψ)}_(m) ⁽²⁾(x_(υj+μ) ^((υ)),y;Δx/υ) without having to evaluate the solutions at any cross range y. These solutions require ten times as much work to compute as {overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx) alone, but this computational expense can be justified by the great increase in range over which the extrapolated solutions are accurate.

We denote these estimates by _(Δx) ⁽⁴⁾∘{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y). They are computed as follows. _(Δx) ⁽⁴⁾∘{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y)=[8195{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx/4)−6561{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx/3) +896{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx/2)−7{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx)]/2520  (9.67) In the modified solution method {overscore (ψ)}_(m) ⁽²⁾is replaced by _(Δx) ⁽⁴⁾∘{overscore (ψ)}_(m) ⁽²⁾. Thus, we use the modified solution method to find the values {overscore (ψ)}_(m)(x_(j+1),y_(k)) and the extrapolates _(Δx) ⁽⁴⁾∘{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k)) for m=1, . . . ,M and k=0,1, . . . ,N from the values {overscore (ψ)}_(m)(x_(j),y_(l)) and from the extrapolates _(Δx) ⁽⁴⁾∘{overscore (ψ)}_(m) ⁽²⁾(x_(j),y_(l)) for m=1, . . . ,M and l=0,1, . . . ,N. These quantities are related by eq. (9.37) after we have replaced {overscore (ψ)}_(m) ⁽²⁾ with _(Δx) ⁽⁴⁾∘{overscore (ψ)}_(m) ⁽²⁾ in it. {{overscore (Ψ)}_(m)(x_(j+1)−0,y_(k))−k_(m) _(Δx) ⁽⁴⁾∘{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y_(k))}_(k=0) ^(N) =ℑ_(N) ⁻¹(exp{−iη²Δx/2k_(m)}ℑ_(N){{overscore (Ψ)}_(m)(x_(j)+0,y_(l))−k_(m) _(Δx) ⁽⁴⁾∘{overscore (ψ)}_(m) ⁽²⁾(x_(j),y_(l))}_(l=0) ^(N)).  (9.68)

We also need the boundary values {overscore (ψ)}_(m) ⁽¹⁾(x_(υj+μ) ^((υ)),±W) at the interpolated range steps to evaluate the boundary conditions that are satisfied by the solutions {overscore (ψ)}_(m) ⁽²⁾(x_(υj+μ) ^((υ)),y;Δx/υ). We compute them by replacing Δx with μΔx/υ in eq.(9.44). Thus, for μ=1, . . . ,υ, k m ⁢ φ ~ m ( 1 ) ⁡ ( x vj + μ ( v ) , ± W ) = N - 1 ( y k = ± W ) ⁢ ( exp ⁢ { - in 2 ⁢ μΔ ⁢   ⁢ x / v2k m } ⁢ N ⁢ { Ψ ~ m ⁡ ( x j + 0 , y l ) - k m ⁢ Δ ⁢   ⁢ x ( 4 ) ∘ φ ~ m ( 2 ) ⁡ ( x j , y l ) } l = 0 N ) . ( 9.69 ) Note that it is not necessary to compute the whole inverse DFT to evaluate the RHS of this equation because we need only the values from the inverse DFT at the endpoints y_(k)=±W.

In medium range applications, the lower order Richardson extrapolates _(Δx) ⁽²⁾∘{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y)=[4{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx/2)−{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx)]/3,  (9.70) _(Δx) ⁽³⁾∘{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y) =[243{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx/3)−128{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx/2)+5{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y;Δx)]/120  (9.71) can be used in place of _(Δx) ⁽⁴⁾∘{overscore (ψ)}_(m) ⁽²⁾(x_(j+1),y) to reduce the computational cost of the procedure. 10. Ramifications.

Many modifications of the invention are possible in light of the above teachings. For example, it can be adapted to waveguides with multiple internal interfaces. In addition, loss can be incorporated in the model. One way to do this is to modify the modal wavenumbers by a perturbative method that is well-known in the art. Furthermore, a variable step size and variable reference wavenumbers can be used in place of the fixed step size and fixed reference wavenumbers in the preferred embodiment.

Additional variations of the invention also are possible in light of the above teachings. Thus, it can be adapted for predicting propagation of seismic waves in the earth's crust, propagation of radio waves in the earth 's atmosphere, and propagation of light waves in dielectric material. Furthermore, computer implementation of the invention is most efficient on parallel computers because many of the operations in the decomposition method can be done in parallel or can be individually vectorized. It is therefore to be understood that the invention may be practiced otherwise than has been specifically described. 

1. A computer implemented method for predicting propagation of a wave in a medium, said medium bordering a top reflective surface and a bottom reflective surface, said method comprising the steps of: (a) selecting a direction in which predictions of said propagation are desired; (b) selecting a pair of transparent boundaries that are parallel to said direction; (c) selecting a plurality of range steps in said direction; (d) calculating a plurality of local modes between said transparent boundaries and at said range steps; (e) calculating a plurality of coupling coefficients between said transparent boundaries and at said range steps; (f) selecting a plurality of modal amplitudes between said transparent boundaries and at said range steps; (g) splitting said modal amplitudes into a plurality of primary components and a plurality of secondary components; (h) selecting a plurality of initial values for said modal amplitudes at a first one of said range steps; (i) using a decomposition method to predict said modal amplitudes between said transparent boundaries and at said range steps from said initial values, said decomposition method including (1) using a coupled-parabolic-equation method with said coupling coefficients to predict boundary values of said primary components at an arbitrary range step from said modal amplitudes and said secondary components at a preceding range step, (2) using a parabolic-equation method to predict said secondary components at said arbitrary range step from said boundary values and from said secondary components at said preceding range step, (3) using said coupled-parabolic-equation method with said coupling coefficients to predict said modal amplitudes at said arbitrary range step from said secondary components at said arbitrary range step and from said modal amplitudes and said secondary components at said preceding range step, j) determining mathematically said propagation of said wave from said local modes and from said modal amplitudes.
 2. A method according to claim 1 in which said wave is a sound wave and in which said medium is an ocean.
 3. A method according to claim 1 in which said wave is a seismic wave and in which said medium is the earth's crust.
 4. A method according to claim 1 in which said wave is a radio wave and in which said medium is the earth's crust.
 5. A method according to claim 1 in which said wave is a light wave and in which said medium is a dielectric material. 